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Abstract 



S^ . We study the relation between the dynamical critical behavior and the kinematics of 

H I stochastic multigrid algorithms. The scale dependence of acceptance rates for nonlocal 

Metropolis updates is analyzed with the help of an approximation formula. A quantitative 
study of the kinematics of multigrid algorithms in several interacting models is performed. 
We find that for a critical model with Hamiltonian 7^(0) absence of critical slowing down 
can only be expected if the expansion of (7^(0 + V')) in terms of the shift ■0 contains no 
relevant term (mass term). The predictions of this rule are verified in a multigrid Monte 
Carlo simulation of the Sine Gordon model in two dimensions. 

Our analysis can serve as a guideline for the development of new algorithms: We pro- 
pose a new multigrid method for nonabelian lattice gauge theory, the time slice blocking. 
For SU{2) gauge fields in two dimensions, critical slowing down is almost completely elim- 
inated by this method, in accordance with the theoretical prediction. The generalization 
of the time slice blocking to SU{2) in four dimensions is investigated analytically and 
by numerical simulations. Compared to two dimensions, the local disorder in the four 
dimensional gauge field leads to kinematical problems. 
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1 Introduction 

Monte Carlo simulations have become an important tool for the study of critical phenomena in 
statistical mechanics |I| and for nonperturbative calculations in Euclidian lattice field theory 
close to the continuum limit 0. 

However, the method has limitations. In the vicinity of a critical point the phenomenon of 
critical slowing down (CSD) is a serious problem: For conventional local algorithms the auto- 
correlation time - that is, roughly speaking, the time needed to generate a new, "statistically 
independent" configuration on a computer - grows rapidly as the system approaches criticality. 
More precisely, the autocorrelation time r behaves like t ~ i^^, where ^ denotes the spatial 
correlation length, and z is the dynamical critical exponent. For conventional local algorithms, 
z ^ 2. Thus, when the critical point is approached, there is a dramatic increase of computer 
time needed to calculate observables to a given accuracy. It is therefore important to find 
Monte Carlo algorithms that have reduced CSD. 

Accelerated algorithms that are still local, such as overrelaxation or the optimized hybrid 
Monte Carlo algorithm, can sometimes reduce the dynamical critical exponent to ^ ~ 1 [^, §|. 

For the complete elimination of CSD in the sense of 2; ~ 0, nonlocal update algorithms are 
needed. Heuristic arguments based on the picture that in simulations with local algorithms 
"information" propagates though the lattice in the form of a random walk lead to a value of 
z = 2. Therefore it is natural to use nonlocal update schemes that perform changes of the 
configuration not only on the scale of the lattice spacing but create fluctuations on all length 
scales, in analogy to the physical fluctuations at a critical point. 

Stochastic cluster algorithms are very successful in reducing CSD in spin models: The first 
cluster algorithm was developed by Swendsen and Wang for discrete spin variables in Ising 
and Potts models |^ . This method could be successfully extended to continuous spin variables 
in 0{N) nonlinear a-models by Wolff ^ and Hasenbusch |^. However, the generalization to 
SU{N) X SU{N) principal chiral models with N > 2, which are models with the same type of 
variables as QCD, seems to be difficult [Q. There are arguments that cluster algorithms will 
not be efficient if the spin variables take values in manifolds other than spheres, products of 
spheres or the quotient of such a space by a discrete group [^ . 

An alternative method to overcome CSD in the simulation of models with continuous vari- 
ables is multigrid Monte Carlo: Multigrid methods are well established tools for the solution of 
discretized partial differential equations [|r^, |rT|. Their generalization to Monte Carlo simula- 



tions was proposed by Parisi |jT2|. Early ideas for nonlocal updating schemes were formulated 



by H. Meyer- Ortmanns ||T^. Goodman and Sokal generalized deterministic multigrid methods 



to multigrid Monte Carlo in the (f) theory in two dimensions fT^ . Mack presented a stochastic 



multigrid approach that was inspired by constructive quantum field theory and renormalization 



group considerations |15 



In the present thesis, every algorithm that performs stochastic updates of predesigned shape 
on a hierarchy of length scales is called multigrid Monte Carlo algorithm. 

There are examples, e.g. the SU{3) x SU{3) principal chiral model, where no fast clus- 
ter algorithms have been found whereas multigrid Monte Carlo algorithms could reduce CSD 
considerably H, |T^ . 

Presently, the only generally applicable method to study the dynamical critical behavior 
of Monte Carlo algorithms for interacting models is numerical experiment. For some models 
experiments show that the dynamical critical exponent z can be substantially reduced by a 



multigrid algorithm [16, |T^, [T^, 19 1. For other models, still z = 2 is found [T^, 20]. Theoretical 
insight into the critical dynamics of multigrid Monte Carlo algorithms is therefore desirable. 

We present a method that can help to judge which algorithm will have a chance to over- 
come CSD in the simulation of a given model before performing the simulation: We study the 
kinematics of multigrid Monte Carlo algorithms. By kinematics we mean the study of the 
scale (block size) dependence of the Metropolis acceptance rates for nonlocal update proposals. 
We do not address the more difficult problem of analytically calculating the dynamical critical 
behavior from the stochastic evolution of the system. However, we are able to find a heuristic 
relation between the kinematical and the critical dynamical behavior of multigrid Monte Carlo 
methods. This relation is based on the fact that sufficiently high acceptance rates are necessary 
to overcome CSD. 

We derive an approximation formula for the block size dependence of acceptance rates for 
nonlocal Metropolis updates. The influence of the coarse-to-fine interpolation kernel (shape 
function) on the kinematics in free field theory, where the formula is exact, is investigated in 
detail. 

The formula is then applied to several interacting models and turns out to be a very good 
approximation. By comparison with free field theory, where CSD is eliminated by a multigrid 
algorithm, we can formulate a necessary criterion for a given multigrid algorithm to eliminate 
CSD: For a critical model with a fundamental Hamiltonian 7i(0) absence of CSD can only be 
expected if the expansion of {7i{(j) + ip)) in terms of the shift ip contains no relevant term (mass 
term). 

Heuristically, the physical content of this criterion can be summarized in a Leitmotiv: 
A piecewise constant update of a nonlocal domain should only have energy costs proportional 
to the surface of the domain, not energy costs proportional to the volume of the domain. 

This seems to be a feature that multigrid algorithms share with cluster algorithms |p . 

As a ffist test of the predictive power of the kinematical analysis we perform a multigrid 
Monte Carlo simulation of the Sine Gordon model in two dimensions. There, the criterion tells 
that we have to expect CSD caused by too small acceptance rates of the updates on large scales. 
This prediction is conffimed by the numerical experiments. In addition it is studied whether 
one can compensate for too small amplitudes of the updates on large blocks by accumulating 
many of these updates randomly. 

After the introduction of the kinematical analysis of multigrid Monte Carlo algorithms in 
the context of spin models, where the multigrid methods have already been developed, we show 
that our method can also be helpful in the design of new multigrid procedures. 

One of the challenges in the development of fast Monte Carlo algorithms is pure lattice 
gauge theory in four dimensions. In particular the nonabelian gauge groups SU{2) and SU{3) 
are important for nonperturbative calculations in the standard model of elementary particles. 

For the dynamical critical behavior of the local heat bath algorithm in SU{3) only very 



crude estimates are available up to now |21], consistent with z ^ 2. The present state-of- 



the-art algorithm for nonabelian gauge fields in four dimensions is overrelaxation [Q . For this 
algorithm, first estimates for z in SU{2) lattice gauge theory gave z = 1.0(1) in physically small 
volumes pH] . Effort was also spent in developing nonlocal algorithms for gauge theories. A fast 



cluster algorithm was found for 3 + 1-dimensional SU{2) gauge theory on a L^ x T lattice at 
finite temperature, but only in the special case T = 1 [^]. A cluster algorithm for ^(1) gauge 
theory in two dimensions is based on the reduction of the gauge theory to a one dimensional 



XY model [25]. However, apart from these special cases, no efficient cluster algorithm for 
continuous gauge groups has been found up to now. 

Therefore possible developments of stochastic multigrid methods for pure gauge fields are of 
particular interest. Multigrid algorithms for U{1) gauge models were introduced and studied in 
two and four dimensions |TB|, ^ . A different but related nonlocal updating scheme for abelian 



lattice gauge theory is the multiscale method ||27| 



We propose a new multigrid algorithm for nonabelian gauge theory and analyze its kine- 
matics. 

To gain experience, we ffist study the case of gauge group SU{2) in two dimensions. We 
introduce a multigrid method for nonabelian gauge theory that treats different time slices 
independently: the time slice blocking algorithm. The theoretical analysis predicts that CSD 
can be eliminated by the time slice blocking. By numerical experiments on systems with lattice 
sizes up to 256^ we check whether this is indeed the case. 

In a second step we generalize the time slice blocking to SU{2) in four dimensions. Compared 
to the two dimensional case we have to face additional difficulties that are caused by the local 
disorder of the gauge fields. Our approximation formula turns out to be very reliable also in this 
case and allows for a prediction of acceptance rates for a large class of nonlocal updates. We 
attempt to estimate the kinematical behavior of the proposed algorithm in the weak coupling 
limit and study whether a reduction of CSD can be expected. In numerical simulations, the 
time slice blocking algorithm in SU{2) in four dimensions is compared with a local heat bath 
algorithm. 

When Parisi proposed the use of multigrid methods for Monte Carlo simulations in 1983, 



he stated |T2|: "On the contrary [to quadratic actions] the application of the multigrid method 
to the gauge field sector seems to be particularly painful but it may be rewarding." 

We attempt to describe where we stand ten years after Parisi's speculation. 

This thesis is organized as follows: In the ffist part (sections |]-0) the kinematical analysis 
of multigrid algorithms is introduced and discussed. In the second part (sections ISHTTl) the 
method of the kinematical analysis is applied to the development of new multigrid algorithms 
for nonabelian gauge fields. 

In section ^ we introduce multigrid Monte Carlo algorithms. Section |^ contains the deriva- 
tion of our approximation formula for acceptance rates. Several coarse-to-fine-interpolation 
kernels are discussed in section §. In section ^ the acceptance rates in free field theory are 
examined in detail. The main idea of the thesis is introduced and discussed in section ^ The 
kinematical analysis for the Sine Gordon, XY, 0"^, 0{N) and CP^~^ models is presented. In 
section |^ the prediction of this analysis for the Sine Gordon model is verified by a multigrid 
Monte Carlo simulation. 

Section ^ starts with the discussion of abelian gauge fields in two dimensions. There are two 
directions of increasing the difficulty of the problem: from abelian gauge group to nonabelian 
gauge group and from two dimensions to four dimensions. Most of the concepts that are needed 
for the treatment of nonabelian gauge fields are introduced in the context of SU{2) lattice gauge 
theory in two dimensions. A new multigrid procedure, the time slice blocking, is introduced. 
In section ^ SU{2) lattice gauge theory in two dimensions is simulated with the time slice 
blocking. 

The new procedure is generalized from SU{2) in two dimensions to SU{2) in four dimensions 
in section |l^. Section ^ contains the results of a multigrid Monte Carlo simulation of SU{2) 



lattice gauge theory in four dimensions. A summary is given in section 112. 



Parts of this thesis have been pubhshed before in refs. [^, ^, |32, |3l|, |32 



2 Multigrid Monte Carlo algorithms 

We consider lattice models with partition functions 

Z= [ U d<f)^ exp{-n{(j))) (2.1) 

on hypercubic c?- dimensional lattices Aq with periodic boundary conditions. The lattice spacing 
is set to one. We use dimensionless spin variables 0^.. An example is single-component 0^ theory, 
defined by the Hamiltonian 

n{<P) = i(0, -A0) + ^E^l + ^E^t^ (2-2) 



2 ^-^^' ■ 4! 



where 



-A0)= E (0X-0,)'. (2.3) 



<x,y> 



The sum in eq. (|2.3| ) is over all nearest neighbor pairs in the lattice[|. 

2.1 Local Metropolis updating 

A standard algorithm to perform Monte Carlo simulations in a model of the type defined above 
is the local Metropolis algorithm: One visits in a regular or random order the sites of the lattice 
and performs the following steps: At site Xo, one proposes a shift 



.^ + s. (2.4) 



The configuration {(px} remains unchanged for x ^ Xq- s is a random number selected according 
to an a priori distribution p{s) which is symmetric with respect to s ^ —s. E.g., one selects 
s with uniform probability from an interval [—£,£]. Then one computes the change of the 
Hamiltonian 

AH = 1-L{<p') - n{<p) . (2.5) 

Finally the proposed shift is accepted with probability min[l, exp(— A7i)]. Then one proceeds 
to the next site. 

The local Metropolis algorithm suffers from CSD when the correlation length in the system 
becomes large: long wavelength fluctuations cannot efficiently be generated by a sequence 
of local operations. It is therefore natural to study nonlocal generalizations of the update 
procedure defined above. 

2.2 Nonlocal multigrid updating 

Consider the fundamental lattice Aq as divided in hypercubic blocks of size l'^, where / denotes 
the coarsening factor. Typical coarsening factors are / = 2 or 3. This defines a block lattice 
Ai. By iterating this procedure one obtains a whole hierarchy of block lattices Aq, Ai, . . . ,Kk 
with increasing lattice spacing (see figure [^). This hierarchy of lattices is called multigrid. 

^The definitions for lattice gauge theory will be introduced in section g 
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Figure 1: Division of the finest lattice in block lattices of Lb x LB-blocks in two dimensions 



Let us denote block lattice points in A^ by x'. Block spins ^x' are defined on block lattices 
Afc. They are averages of the fundamental field (p^ over blocks of side length Lb = l^'- 



^x' 



Xdx' 



(2.6) 



The sum is over all points x in the block x' . The L^-dependent factor in front of the average 
comes from the fact that the corresponding dimensionful block spins are measured in units of 
the block lattice spacing: A scalar field 0(x) in d dimensions has canonical dimension {2 — d)/2. 



Thus (f){x) 



,(2- 



^^^'^(px, where a denotes the fundamental lattice spacing. Now measure the 



dimensionful block spin $(x') in units of the block lattice spacing a': $(x') = a'^'^^'^^^'^^^' , with 
a' = uLb- If we average in a natural way $(x') = L^ J2xi^x' (l^{^) a-^id return to dimensionless 
variables, we obtain eq. (W^ )- 

A nonlocal change of the configuration (p consists of a shift 



4>x ^ 4>x + S^Jx 



(2.7) 



s is a real parameter, and the "coarse-to-fine interpolation kernel" (or shape function) ipj. 
determines the shape of the nonlocal change. ip is normalized according to 



Lb' E ^^ 

x£x' 



L 



(2-d)/2 
B 



Sx'x' 



(2. 



s for x' = X' 



Note that by the nonlocal change ( |2.7] ), the block spin is moved as ^x' -^ ^x' 
and remains unchanged on the other blocks. The simplest choice of the kernel i(j that obeys the 
constraint ( |2.8D is a piecewise constant kernel: ipx = Lb , if x G x'^, and otherwise. Other 
kernels are smooth and thus avoid large energy costs from the block boundaries. A systematic 
study of different kernels will be given in section |^. Two one dimensional examples for ip are 
given in figure |^. 
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Figure 2: Coarse-to-fine interpolation kernels ip in one dimension: Top: piecewise constant 
kernel, bottom: piecewise linear kernel 



2.3 The average Metropolis acceptance rate Q(s) 

The s-dependent Metropolis acceptance rate for such proposals is given by 

n{s) = (min[l, exp{-An)]) . (2.9) 

Here, ((.)) denotes the expectation value in the system defined by eq. (|2.1|). Furthermore, 

AH = n{(l) + sij) - n{(f)) . (2.10) 

Q{s) can be interpreted as the acceptance rate for shifting block spins by an amount of s, 
averaged over a sequence of configurations generated by a Monte Carlo simulation. As a static 
observable, Q{s) is independent of the algorithm that we use to compute it. Q{s) is a useful 
quantity when one wants to know whether updates with increasing nonlocality (i.e. increasing 
block size Lb) can be performed in an efficient way. Of course, different choices of the kernel 
ijj result in different acceptance rates. 

In actual Monte Carlo simulations, s is not fixed. In the same way as in the local Metropolis 
algorithm, s is a random number distributed according to some a priori probability density. If 
we choose s to be uniformly distributed on the interval [—£,£:], the integrated acceptance rate 
Pace (as customarily measured in Monte Carlo simulations) is obtained by averaging Q{s) as 
follows: 

Pacc{e) = ^£jsn{s) . (2.11) 

It turns out to be a good rule to adjust the maximum Metropolis step size e such that 



P. 



0.5. 



2.4 Unigrid versus recursive multigrid 

We consider every algorithm that updates stochastic variables on a hierarchy of length scales 
as multigrid Monte Carlo algorithm. However, there are two different classes of multigrid algo- 
rithms: multigrid algorithms in a unigrid implementation and "recursive" multigrid algorithms. 



In the unigrid formulation one considers nonlocal updates of the form (|2.7| ). Updates on the 
various layers of the multigrid are formulated on the level of the finest lattice Aq. There is no 
explicit reference to block spin variables $ defined on coarser layers A^ with k > 0. In addition, 
unigrid also refers to a computational scheme: Nonlocal updates are performed directly in 
terms of the variables on the finest grid Aq in practical simulations. An example for a unigrid 
implementation is explained for two dimensional nonabelian gauge theory in section p. 

In contrast, the recursive multigrid formulation consists of recursively calculating conditional 
Hamiltonians that depend on the block spin variables $ on coarser layers A^. This formulation 
is possible if the conditional Hamiltonians are of the same type or similar to the Hamiltonian 
on the finest lattice. Then, the conditional probabilities used for the updating on the fc-th layer 
can be computed without always going back to the finest level Aq. Therefore, an recursive 
multigrid implementation reduces the computational work on the coarser layers (see the work 
estimates below). At least in free field theory, a recursive multigrid implementation with smooth 
interpolation is possible using 9-point prolongation kernels (a special case of piecewise linear 
kernels as shown in figure |^) in two dimensions and generalizations thereof in higher dimensions 
P^l , [5^ . Generally, a recursive multigrid implementation for interacting models is only feasible 
in special cases with piecewise constant kernels. Examples for such implementations are given 
for the Sine Gordon model in section |^ and in the four dimensional nonabelian gauge theory in 



section |TT] below. 

An algorithm formulated in the recursive multigrid style can always be translated to the 
unigrid language (that is the way we are going to use the unigrid formulation for the analysis 
of multigrid algorithms). The reverse is not true, since not all nonlocal changes of the field 
on the finest lattice can be interpreted as updates of a single block spin variable of a recursive 
multigrid. As an example, one can use stochastically overlapping blocks in the unigrid style by 
translating the fields by a randomly chosen distance [|l^] . 

If we formulate our kinematical analysis in the unigrid language we nevertheless can include 
all algorithms formulated in the recursive multigrid style. 

2.5 Multigrid cycles 

The sequence of sweeps through the different layers A^ of the multigrid is organized in a periodic 
scheme called cycle [|l5]. The most important cycles are illustrated in figure |^. The simplest 
scheme is the V-cycle: The sequence of layers visited in turn is Aq, Ai, . . . A^, Ax_i . . . Ai. 
More general cycles are characterized by the cycle control parameter 7. The rule is that from 
an intermediate layer A^ one proceeds 7 times to the next coarser layer A^+i before going back 
to the next finer layer Afc_i. A cycle control parameter 7 > 1 samples coarser layers more 
often than finer layers. With 7 = 1 we obtain the V-cycle. 7 = 2 yields the W-cycle that is 
frequently used with piecewise constant kernels. 



The computational work estimates for the different cycles are as follows [Q: The work 
for a recursive multigrid cycle is ~ L'^ if 7 < Z'^, where L denotes the lattice size and / is the 
coarsening factor. The work for a unigrid cycle is ~ L'^logL if 7 = 1, and ~ L'^+'^°^i'y if 7 > 1. 

8 





Figure 3: V-cycle ('-/ = 1) and W-cycle ('-/ = 2) 



If one wants the computational work in the unigrid style not to exceed (up to a logarithm) 
an amount proportional to the volume L'^ of the lattice, one is restricted to use a V-cycle. 
Simulations with 7 > 1 (e.g. a W-cycle) can only be performed in the recursive multigrid style. 



3 An approximation formula for Q(s) 

In this section we shall derive an approximate formula for the quantity Q{s) defined in ( p.9|) . 
We can write Q{s) as 

n{s) = [du min(l, e-") / — e'^^" (e^^^^) . (3.1) 

Let us assume that the probability distribution of A7Y is approximately Gaussian. We param- 
eterize this distribution as follows: 

dprob(A7^) oc dAH exp( — ^{AH - h^f) , (3.2) 

with hi = (AH) and /12 = K(^^^) - (^^)^)- Then 

(exp(ipA7-^)) ^ exp{ihip — h2P^) . (3.3) 

We now show that from the assumption that the probability distribution of A7i is Gaussian 
we can derive that hi = h2- The starting point is the identity 

1 = ^/ n d<P,exp{-nm. (3.4) 

xGAo 

If we use the translational invariance of the measure OxeAo ^^^ — YixeAo difpx + sipx) we get 

^ = \ I X{d4>x exp{-n{(j) + #)) = ^ f Y[d<f>. exp{-n{<j) + #) + H{(j)) - H{<j))) . (3.5) 
With the definition ATi = 7i(0 + sip) — 7i(0) we obtain 

1 = 4/ n #- exp{-AH)exp{-n{<P)) = {exp{-An)) . (3.6) 

This identity is independent of the form of the Hamiltonian. If we now set p = i in eq. ( |3.3| ) 
and compare this with eq. ( p.6|) , we are led to hi = h2. Therefore eq. ( p.3|) reduces to 

{expiipATi)) :^ exp[(ip — p^)/ii] , (3.7) 

and the integrations in eq. (p.l|) can be performed exactly since there are only Gaussian integrals 
involved. The result is 

n{s) ^ eiic{\^i) . (3.8) 

with erfc(x) = 2/y^/^t/t exp(— t^). (For an analogous result in the context of hybrid Monte 
Carlo see 11.) 

For free massless field theory with Hamiltonian 7i(0) = ^(0, —Acj)), we get hi = \a s^ with 
a = {ip, —Aip), and our approximation formula becomes exact 

fi(.)=erfc(y||.|). (3.9) 



Eq. (|3.9|) can be checked directly by using (exp(ipA7i)) = exp[{ip — p'^)hi] in eq. ( p.l|) . This 
relation is exact in free field theory. 
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4 Coarse-to-fine interpolation 

In this section we shall discuss several choices of the coarse-to-fine interpolation kernels. In 
order to have a "fair" comparison, all kernels ip will be normalized according to eq. 



4.1 Optimal interpolation kernels 

In free massless field theory, the quantity a = {ip, —Aip) characterizes the decrease of the 



acceptance rate n{s) as given in (3^) with increasing shift s. Therefore it is natural to minimize 
a in order to maximize Q{s) for fixed s. 

The optimal kernel ip^^"-'^^ from the point of view of acceptance rates can be defined as 
follows: minimize the quadratic form 

a = (V',-AV') (4.1) 

under the constraints that the average of?/' over the "central block" x'^ is given by L^^ , and 
its average over blocks x' ^ x'^ vanishes: 

Lb" E V^- = L^b''^'^^^',.'^ for all x' e Ku • (4.2) 

This variational problem can be solved with the help of Fourier methods. The result is 

^e.act ^ Lg+'^)/2^^^^, , (4.3) 

where Ax,x'^ denotes the Gaw§dzki-Kupiainen kernel (see, e.g. [0). The use of this kernel 
leads to a complete decoupling of the different layers of the multigrid. This way of interpolating 
from a coarser block lattice A^ to the fine lattice Aq is well known in rigorous renormalization 
group theory |^. It is interesting that considerations about optimizing acceptance rates in a 



stochastic multigrid procedure lead to the same choice of the interpolation kernel. 

Because ip'^^"-'^^ is nonvanishing on the whole lattice, it is not convenient for numerical pur- 
poses. For an attempt to change the block spin $3./^ on block x'^ one has to calculate contri- 
butions to the change of the Hamiltonian from all lattice points. Therefore the computational 
work for a single update is proportional to the volume. 

4.2 Truncation of the support of the optimal kernel 

4.2.1 Kernels with support on a block and its nearest neighbors 

We define a "truncated kernel" ■j/)*™"'= by restricting the support of ip on the block x'^ and its 
nearest neighbor blocks y'^ 

^trunc ^ g jf a; ^ Xq or X ^ y'^, where y'on.n. x'^ . (4.4) 



In other words, the Laplacian in eq. ([4.1|) is replaced by a Laplacian A^ with Dirichlet boundary 
conditions on the boundary of the support of ip. We again minimize a = {ip, —Aoip) under the 
2d + 1 constraints that the average of ip over the blocks x'^ and its nearest neighbor blocks is 
given. This minimization can be performed numerically by a relaxation procedure. In order to 
maintain the normalization condition, one always updates simultaneously two spins residing in 
the same block, keeping their sum fixed. The '?/'*™"''^-kernels were used in a multigrid simulation 



of the model in four dimensions [^| 
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4.2.2 Kernels with support on a single block 

From a practical point of view, it is convenient to use kernels that have support on a single 
block Xq, i.e. 

^lj, = Q lix^x',. (4.5) 

We define a kernel -i/;™-™ with this property by minimizing a = {ip, —Ao,x'„ip) under the con- 
straint that the average of ip over the block x'^ is given. The Laplacian with Dirichlet boundary 
conditions on the boundary of x'^ is defined as follows: 



{^D,x',(p)x 



-2d<j),+ Y. 



'ry 
yn.n.x 



for X ^ x'g . (4.6) 



^mm ^-.g^j^ i^g calculated using an orthonormal set of eigenf unctions of Ad,x'^- 

4.3 Other kernels with support on the block 

We shall now discuss other kernels with support on the block that are frequently used in the 
literature. 

4.3.1 Piecewise constant interpolation 

Piecewise constant interpolation kernels are defined by 

f T {2-d)/2 r „ ^ „/ 

^--*= ^B lorxfcx ^^^^^ 

This kernel has the advantage that for many models the conditional Hamiltonians used for 
updating on coarse lattices are of the same type or similar to the Hamiltonian on the finest 
lattice. This means that the conditional probabilities used for the updating on the fc-th layer 
can be computed without always going back to the finest level Aq. Therefore, an recursive 
multigrid implementation with cycle control parameters 7 > 1, e.g. a W-cycle, can be used. 

Piecewise constant kernels in a recursive multigrid implementation with a W-cycle will be 
used in the Sine Gordon model in section |^ and in the four dimensional nonabelian gauge theory 



in section 11 



4.3.2 Piecewise linear interpolation 

We consider the block 

x', = {x\x^ G{l,2,3,...,LB},/x = l,...,rf}. (4.8) 

The kernels for other blocks are simply obtained by translation. For Lb even, ip'-^^'^'^^ is given 
by 



d 



A^ is a normalization constant. Piecewise linear kernels are going to be used in the simulation 
of two dimensional nonabelian gauge theory in section |n|. There, smooth nonlocal heat bath 
updates can be performed with linear interpolation. 
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Table 1: Results for a = {ip, —Aip) in 2 dimensions, 512^ lattice 



kernel 


Lb=2 


Lb=4: 


Lb=8 


Lb=16 


Lb=32 


Lb=64 


Lb=128 


Lb=256 


exact 


6.899 


8.902 


9.705 


9.941 


10.00 


10.02 


10.18 


13.11 


trunc 


7.000 


9.405 


10.73 


11.38 


11.69 


11.84 


11.92 


- 


min 


8.000 


13.24 


18.48 


22.58 


25.23 


26.76 


27.59 


28.02 


sine 


8.000 


13.62 


19.34 


23.78 


26.62 


28.25 


29.13 


29.58 


linear 


8.000 


15.80 


24.58 


31.84 


36.68 


39.51 


41.05 


41.84 


const 


8.000 


16.00 


32.00 


64.00 


128.0 


256.0 


512.0 


1024 



4.3.3 Ground state projection 

Ground state projection kernels are defined as follows: ijj^^"-'^ is the eigenfunction corresponding 
to the lowest eigenvalue of the negative Laplacian with Dirichlet boundary conditions — A^i j.^ : 



^: 



sine 

X 



^fY[ Si 



TT 



sm 



M=i 



Lb + 1 







x^) for X E x'g 



for X ^ Xq . 



(4.10) 



Again, J\f denotes a normalization constant. Note that this kernel is different from t^'"*". A 
generalization of this kernel was introduced for scalar fields in the background of nonabelian 
gauge fields in ref. 



4.4 Results for a = {i/j, —Ai/j) in two and four dimensions 

The results for the quantities a = {ip, —Aip) for different kernels in two dimensions are presented 
in table |I|. We used a 512^ lattice (7^^^'"'=* depends on the lattice size). The different kernels are 
ordered according to increasing value of a. 

are close together. This shows that the truncation of the 



The values of q/^^'"'^* and 



a 



trunc 



support of ijj to the block and its nearest neighbor blocks is a good approximation to ijj^^'^'^ (in 
the sense of acceptance rates). The value of a^^'^^ for Lb = 256 is remarkably higher than on 
smaller blocks. This is a finite size effect because the block lattice consists only of 2^ points. 



Since the nearest neighbors overlap on a 2^ lattice, no result for 



a 



trunc 



is quoted for Lb = 256. 



are 

mm T^ 



The values of a for the smooth kernels with support on the block ^™«"^^'5«»^e and ■j^'*"-^'^^ 
of the same magnitude. We can see that %p'^^^'^ is almost as good as the optimal ip = ijj 
contrast to all smooth kernels, a™"^* grows linear in Lb- 

The results for different kernels in four dimensions are presented in table 0. Here we used 
a 64^ lattice. In principle, the a's behave as in two dimensions. The values of q/'*"^^"'" for small 
blocks are higher than a*^""**. The pyramids of the piecewise linear kernels have a lot of edges 
in four dimensions which lead to high costs in the kinetic energy. 
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Table 2: Results for a = {ip, —^ip) ^^ 4 dimensions, 64'^ lattice 



kernel 


Lb=2 


Lb=4: 


Lb=8 


Lb=16 


Lb=32 


exact 


14.48 


20.38 


23.48 


24.71 


30.62 


trunc 


14.67 


21.61 


26.54 


29.26 


- 


min 


16.00 


27.72 


41.56 


54.18 


63.33 


sine 


16.00 


30.37 


48.46 


64.85 


76.44 


linear 


16.00 


39.02 


70.78 


101.0 


122.9 


const 


16.00 


32.00 


64.00 


128.0 


256.0 



The L^-dependence of the a's in c? dimensions is 



a 
a 



is»l 



2dLB for piecewise constant kernels . 
const for smooth kernels . 



TT 



As an example, the expression for a''*"'^ in d dimensions is 

^sine ^ Lf\LB + l)'^2'^+2t/sin^'^+2 

For large block sizes we find 



2(L 



B 



sm 



-2d 



n 



w 



a 



Lb»1 



d 



n 



2d+2 



23d 



const 



(4.11) 



(4.12) 



(4.13) 



From table |l| we observe that in two dimensions a becomes almost independent of Lb for the 
smooth kernels if the block size is larger than 16. In four dimensions (table ^, we find a{LB) ~ 
const only for 0;^^""^*. The other a's for the smooth kernels have not become independent of Lb 
for the block sizes studied. 
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5 Acceptance rates in free field theory 

In this section we discuss the scale dependence of Metropohs acceptance rates of muhigrid 
Monte Carlo algorithms in free field theory. The kinematical behavior can be related to the 
dynamical critical behavior of the algorithms. This will be the starting point for the theoretical 
analysis of interacting models in section p. In free field theory CSD is eliminated by a multigrid 
algorithm. 

5.1 Massless free field theory 

Recall the exact result (|3.9|) 

n{s) =erfc(^^||s|j 

in massless free field theory. Q{s) is only a function of the product as'^. If we increase the block 
size Lb, the quantity a increases as a function a^Ls) of the block size. In order to keep Q{s) 



fixed when Lb is increased, we have to rescale the changes s like l/Ja{LB)- As a consequence, 
to maintain a constant acceptance rate in massless free field theory, s has to be scaled down 
like 1/\/Lb for piecewise constant kernels, whereas for smooth kernels the acceptance rates for 
large Lb do not depend on the block size. 

Note that this behavior of the acceptance rates for large Lb is not yet reached in four dimen- 
sions for the block sizes studied (except for ^^^'^ct-j g^g ^^^^ ^j^g discussion of the Metropolis 
step size below. 

Let us relate the scale dependence of the Metropolis amplitude s with the dynamical critical 
behavior. For smooth kernels we have s ~ const for large Lb- This indicates that the fluctua- 
tions generated by the algorithm are of the same size on all length scales. Multigrid simulations 
of free massless field theory with smooth interpolation and a V-cycle yield z = I^Q], Q . 

For piecewise constant kernels we have to scale down the Metropolis amplitude s like 
s ~ l/y/Ls. In simulations with piecewise constant interpolation and a V-cycle, z = 1 
was found |^, ^ ^ . At least for free field theory, this disadvantage of the piecewise constant 
kernels can be compensated for by using a W-cycle instead of a V-cycle, resulting in z = 



Smooth kernels in a unigrid implementation can be used only in V-cycle algorithms. An 
exception are 9-point prolongation kernels in two dimensions and generalizations thereof in 
higher dimensions. They can also be used in a recursive multigrid implementation with a 
W-cycle, at least in free field theory (cf. section |^). 

5.2 Scale dependence of the Metropolis step size 

We now illustrate what this rescaling of s means for the Metropolis step size e in an actual 
multigrid Monte Carlo simulation. Let us inspect the integrated acceptance probability defined 
in eq. ( |2.11| ) . If we insert the exact result ( |3.9| ) for massless free field theory 



Pacc{e) = j-J'^dsQ{s) = j-J'^dseTic(^^\s\] , (5.1) 
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Figure 4: Top: Metropolis step sizes s{Lb) for massless free field theory in two 
dimensions, 512"^ -lattice. Bottom: s{Lb) for massless free field theory in four 
dimensions, 64^-iattice. s{Lb) is choosen in such a way that always Pace = 0.5 
holds. Symbols: full circles: ip<^^°-'^* ^ full triangles: ip^^'^"-'^^ empty circles: -^Z;™*", 
empty triangles: ■j/'**"''^, full squares: ■?^'*"'="''^ empty squares: ip'^°^^^. Lines are only 
drawn to guide the eye. 
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we get 



P., 



erfc ( \ —e 



(5.2) 



Pace is only a function of the product ae^. In order to keep Pace fixed (to, e.g. 50 percent) we 



have to rescale s^Lb) hke l/-i/a(Ls), exactly in the same way as we had to rescale s to keep 
VL{s) fixed. This Ls-dependence is plotted in figure 0. 



5.3 Intuitive random walk picture of the W-cycle 

In free field theory it is possible to overcome the the disadvantage s^Lb) ~ 1/^/Lb of the 
piecewise constant kernels compared to elLs) ~ const for smooth kernels by using a W-cycle 
instead of a V-cycle. This can be made plausible by a simple random walk argument. 

A constant accumulated step size on all length scales can be achieved in the following 
way: For step sizes scaling down like s^Lb) ~ 1/^/Lb and a coarsening by a factor of two, the 
Metropolis step size on a next coarser lattice is too small by a factor of l/-\/2- If we assume that 
subsequent update steps within a multigrid cycle are statistically independent and that these 
steps accumulate in a random-walk like way, we can expect to compensate for this decrease 
by increasing the number of updates on the next coarser grid by a factor of two. This can be 
achieved by a higher cycle with cycle control parameter 7 = 2. Since higher cycles are defined 
recursively by going from an intermediate block lattice 7 times to the next coarser lattice before 
going back to the next finer lattice, 7 times more updates on each coarser lattice are performed. 

In free field theory, the assumptions of this argument seem to be correct, and the random 
walk argumentation can explain that for piecewise constant kernels and a W-cycle CSD is 
eliminated. 



5.4 Massive free field theory 

We now discuss massive free field theory with Hamiltonian H(0) = ^(0, [—A + m^]0). We find 
hi = |«m s"^, with am given by 



«m = (^, [-A + m^jip) = a + w? ^ ipl. 

xeAo 



(5.3) 



Therefore the exact result is Vl{s) = erfc (Jam/8 \s\]. A term J2x'4'l scales ~ L^ in arbitrary 
dimensions. For piecewise constant kernels 



E 

xgAo 



const 



For '^''^"^'^-kernels we find 



xeAo 



and for large block sizes 



d „i^'id 



2U|+'^(Lij + l)''sin 



7"2 



IT 



2{Lb + 1) 



sm 



-2d 



TT 



^B 



E 

xgAo 






(5.4) 



(5.5) 



(5.6) 
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If the block size Lb is smaller than the correlation length ^ = 1/m, hi is still dominated by the 
kinetic term s^{il), —Aip), and the discussion is the same as in the massless case. 

As soon as the block size Lb becomes larger than ^, hi is dominated by the "mass" term 
s'^m'^Y^x'^x ~ s'^L'bi ^^^d s has to be rescaled like s ~ I/Lb in order to maintain constant 
acceptance rates. Of course this is a dramatic decrease for large block sizes compared to 
s ~ const in the massless case (using smooth kernels). Block spins on large blocks are essentially 
"frozen" . However this is not a problem for the performance of the algorithm in massive free field 
theory: The effective probability distribution for the block spins $ is given by exp(— Tieffl"^)), 
where 7ieff($) denotes the effective Hamiltonian in the sense of the block spin renormalization 
group P3[ . The physical fluctuations of the block spins are dictated by an effective mass term 



mis E ^l' with ml^ ~ m'Ll ■ (5.7) 



Thus, the algorithmic fluctuations (described by the mass term m? X]^ V^l ~ m?Ly) and the 
physical fluctuations (described by the effective mass ~ ni?L\) behave similarly, and the multi- 
grid algorithm is able to create fluctuations just of the size that is needed by the physics of the 
model. Moreover there is no need to do updates at length scales larger than ^ in order to beat 
CSD. 

In this sense, the discussed algorithmic mass term m? X^x '4^1 is well behaved for free field 
theory, since it decreases with the physical mass in the vicinity of the critical point. Stated 
differently, in massive free field theory the algorithmic mass term scales with the physical mass. 
As we shall see in section ^ for interacting models close to criticality, a different scenario is 
possible. There, it can happen that an algorithmic mass term ~ Y.x V'x persists, whereas the 
physical mass vanishes. If this happens, the multigrid algorithm is not able to produce the 
large critical fluctuations required by the physics, and we can not expect that CSD will be 
eliminated. 

5.5 Behavior of a term Ex ^^ 

The Ls-dependence of a term Y,x i^t will also be needed in the study of the (p'^ theory in section |] 
below. In d dimensions such a term scales ~ L^~ : For piecewise constant kernels 

E {i^r^'Y = L%-' , (5.8) 

xeAo 

whereas using '?/''^*"'^-kernels we find 

,4 



^ (^^--*j =6''L%+^''{LB + lYsm' 

xGAo 



IT 



2iLB + l] 
In the limit of large block sizes this term behaves like 



sin-^-^ ( -^-^ ) . (5.9) 



igAo 




s(*n ,75, Ui^"- (5,10) 
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5.6 Degree of relevance 

In order to summarize the different large-L^-behavior of local operators in the kernel ip, let us 
introduce the degree of relevance in the sense of the perturbative renormalization group: The 
(superficial) degree r of relevance of a local operator in ip which is a polynomial of m scalar 
fields with n derivatives is defined by r = d + m{2 — d) /2 — n. This definition is valid for smooth 
kernels. For large Lb, an operator with degree of relevance r behaves like L^. An operator is 
called relevant if r > 0. As we have seen in the examples above, a mass term has r = 2, and a 
i/j^-teim has r = 4 — d. A kinetic term a = {ip, —Aip) has r = for smooth kernels. 

The only difference for piecewise constant kernels is that a kinetic term behaves like a = 

(^, -A^) oc Lb- 
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6 Acceptance rates for interacting models 



In this section, we shall apply formula ( |3.8| ) in the discussion of multigrid procedures for various 
spin models in two dimensions: the Sine Gordon model, the XY model, the single- component 0'^ 
theory, the 0{N) nonlinear cr-model and the CP^~^ model. The scale dependence of acceptance 
rates for interacting models will be compared with the kinematical behavior in free field theory, 
where CSD is known to be eliminated by a multigrid algorithm. A necessary criterion for the 
successful acceleration of a simulation by a multigrid algorithm is formulated. 

6.1 The two dimensional Sine Gordon model 

The two dimensional Sine Gordon model is defined by the Hamiltonian 

n<P) = ^(0, -A0) - C E cos 27r<j), . (6.1) 

From the point of view of statistical mechanics, this system can be considered as a two 
dimensional surface in a periodic potential. The model exhibits a Kosterlitz-Thouless phase 
transition at PdC)- In the limit of vanishing fugacity (, (3c takes the value 2/7r = 0.6366 .... 
For (3 > (3c the model is in the massless (rough) phase. The fluctuations of the surface are given 
by the surface thickness 

^^ = ((0.-0)^ , (6.2) 



where denotes the average of the field over the lattice. In the massless (rough) phase, the 
surface thickness a^ scales with logL |H3]. In this phase, the cosine-term of the Hamiltonian is 



irrelevant, and the flow of the effective Hamiltonian (in the sense of the block spin renormaliza- 
tion group) converges to that of a massless free field theory: the long distance behavior of the 
theory is that of a Gaussian model. Since multigrid algorithms have proven to be efficient in 
generating long wavelength Gaussian modes, one might naively conclude that multigrid should 
be the right method to fight CSD in the simulation of the Sine Gordon model in the massless 
phase. But this is not so: 

The energy change of a nonlocal update 0x ^ 0x + sip^ is 



s^ , , ... s 



AH = n{<j> + si(j) - H(0) = — (V', -A^) + -(^, -A0) 
— CE {'^os(27r0a;) [cos{27is4'x) ~ 1] ~ sin(27r0a;) sin(27rs'0a;)} . (6.3) 

X 

The Hamiltonian is globally invariant under (p^ -^ —<px- Therefore, (0^;) = and (sin27r0a;) = 
on finite lattices. Using this, we find for the average energy change hi = (A7i) the expression 

hi = ^s' + CCJ2[^ - cos(27r#.)] , (6.4) 

with C = (cos(27r0^)). Recall that hi is the quantity that determines the acceptance rates Q{s) 
by eq. Q: 

n(s) ^ erfc(W/ii) . 



The essential point is that the second term in ( |6.4D is proportional to the block volume Lg 
for piecewise constant and for smooth kernels (cf. the discussion in section ^). This can be 
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Figure 5: Q{s) for piecewise constant kernels in the two dimensional Sine Gor- 
don model, (3 = 1.0, ( = 1.0. From top to bottom: Lb = 4, 8, 16, 32 on a 
16^, 32^, 64^, 128^ lattice, respectively. Points with error bars: Monte Carlo re- 
sults, lines: analytical results 
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checked for small s by expanding the cosine term in s. One therefore has to face a dramatic 
decrease of acceptance when the blocks become large, even for small fugacity (. A constant 
acceptance rate is achieved only when the proposed update steps are scaled down like I/Lb- 
It is therefore unlikely that any multigrid algorithm - based on nonlocal updates of the type 
discussed here - will be successful for this model. This prediction will be verified by a multigrid 
Monte Carlo simulation of the two dimensional Sine Gordon model in section |^. 

We demonstrate the validity of formula ( |3.8| ) (using a Monte Carlo estimate for C) by 
comparing with Monte Carlo results at (3 = 1.0, C = 1.0. This point is in the massless phase. 
There, the physical scale of the system is set by the linear size L of the lattice. In figure |^ we 
show both the numerical and analytical results for Q{s) for Lb = 4, 8, 16, 32 on lattices of size 
16^ 32^, 642andl282, respectively. 

We tested the precision of our approximation formula for piecewise constant kernels only. 
However, we have no doubts that the quality of the approximation is also very good for other 
interpolation kernels ip. 

6.2 The two dimensional XY model 

We now discuss the two dimensional XY model, defined by the partition function 

Z= [Ude, exp(/3 Y. cos(0. -9y)). (6.5) 

X <x,y> 

The sum in the exponent is over all unordered pairs of nearest neighbors in the lattice. As 
the Sine Gordon model, the XY model has a massless (spin wave) phase for (3 > jSc, and a 
massive (vortex) phase for (3 < (3c- The best available estimate for the critical coupling is 
/?, = 1.1197(5) [H. 

Nonlocal updates are defined by 

Ox ^ Or, + sipx , (6.6) 



with ^p obeying again the normalization condition ( ^.8| ). To define a (linear) block spin, we 
rewrite the partition function ( |6.5|) in terms of 2-component unit vector spin variables Sx'- 

Z= JW {dlsx 6isl - 1)) exp(/3 J2 ^- ■ ^y) (6-7) 

X <x,y> 

The block spins Sx' are then defined as block averages of the unit vectors Sx- 

Note that the proposal (|6.6|) changes the block spin by an amount ~ s only when the spins 
inside the block are sufficiently aligned. This will be the case in the spin wave phase for large 
enough f3. For smaller (3, the correct (or "fair") normalization of the kernels ■?/' is a subtle point. 
We believe, however, that our argument is not affected by this in a qualitative way. 

Let us inspect the energy change of the nonlocal update (^7^ 



An = -I3 J2 {cos(^x -Oy + si^x - ^y)) - cos(^, - Oy)} 



<x,y> 



-P XI {cos(6'a,. - 6y) [cos{s{i)x - ^y)) " 1] - sin(6'^ - Oy) sin(s(^^ - iljy))} . (6.8) 

<x,y> 
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Since the Hamiltonian is globally invariant under 9^ -^ —dxi the expectation value {sm{9x — 0y)) 
vanishes on finite lattices. Therefore, the relevant quantity hi is 

hi = l3E Y.[l- cos(s(7/;, - ^y))] , (6.9) 

<x,y> 

with E = {cos{9x — 9y)), x and y nearest neighbors. For piecewise constant kernels, hi is 
proportional to Lb- For smooth kernels hi will become independent of Lb for large enough 
blocks. For small s, 

hi ^ ^s^(3E Y. {ijj, - ^yf = y^(3Ea . (6.10) 

<x,y> 

As above, a = {ip,—Aip). This quantity becomes nearly independent of Lb already for Lb 
larger than 16 (cf. section ^). 

From the point of view of acceptance rates the XY model therefore behaves like massless 
free field theory. A dynamical critical exponent z consistent with zero was observed in the 
massless phase [^. In this phase, the physically dominant excitations are spin waves. From 
this comparison we conclude that the analysis of the kinematics of multigrid algorithms can 
tell us whether spin wave excitations can be accelerated by the algorithm. 

In the massive phase the physically important excitations are vortices, which are topo- 
logically nontrivial objects. In this phase multigrid Monte Carlo simulations with piecewise 
constant kernels yielded z ^ lA [|^. This is an example for the fact that good acceptance 
rates are not sufficient to overcome CSD. The kinematical analysis is not sensitive to topological 
problems. 

We again checked the accuracy of formula (|3.8| ) by comparing with Monte Carlo results 
at (3 = 1.2 (which is in the massless phase). The results are displayed in figure ^ The only 
numerical input for the analytical formula was the link expectation value E. 

6.3 The 0^ theory in d dimensions 

Let us now turn to single- component rf- dimensional 0^ theory {d = 2,3,4), defined by the 

Hamiltonian 

fji-^ \ 

H(0) = i(0,-A0) + -^E'/'x. + 7TE'/>'- (6-11) 

■^ X ^- X 

A nonlocal update 0x — ^ 0x + sip^ leads to the energy change 

AH = n{(f) + s^) - n{(l)) = sm^j, -A^) + s{^, -A0) 

+ ^ E [^V' + 25^.0.] + ^ E [^V' + 43^^^^. + Qs'^'jl + As^Px<|)l\ ■ (6-12) 

^ X ^- X 

Since ?-^(— 0) = 7-^(0), the expectation values (0i^) and (02.-) vanish on finite lattices. Therefore 
we find for hi = (ATi) 



hi = s^ {la + 



2 4 



E^4 + ''lf^^'' ^^-^^^ 



where P = {(pl). Recall that Y.xi'l increases with L^, independent of d, whereas Y.xi't scales 
like L^~ , for smooth and for piecewise constant kernels (cf. the discussion of the different 
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Figure 6: Q(s) for piecewise constant kernels in the two dimensional XY model, 
P = 1.2. From top to bottom: Lb = 4, 8, 16 on a 16^, 32^, 64^ lattice, respectively. 
Points with error bars: Monte Carlo results, lines: analytical results 
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Figure 7: Q{s) for piecewise constant kernels in the two dimensional 0^ theory, 
ml = -0.56, \o = 2.4. From top to bottom: Lb = 4, 8, 16 on a 16^, 32^, 64^ 
lattice, respectively. Points with error bars: Monte Carlo results, lines: analytical 
results 
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kernels in section ^. We conclude that also in this model we have to face rapidly decreasing 
acceptance rates when the blocks become large. As in the case of the Sine Gordon model, s has 
to be rescaled like 1/Lb in order to maintain constant acceptance rates. 

Therefore there is little hope that multigrid algorithms of the type considered here can 
overcome CSD in the single-component (p^ model. In numerical simulations of two dimensional 
0^ theory, a dynamical critical behavior is found that is consistent with z ^ 2 for piecewise 
constant and for smooth kernels [0, 0, ^. In four dimensions, there is no definite estimate 
for z 0. 

Figure ^ shows a comparison of Monte Carlo results for two dimensional 0^ theory with the 
theoretical prediction based on the numerical evaluation of P. The simulations were performed 



in the symmetric phase at m^ 



-0.56 and Xq = 2.4. The correlation length at this point is 



The quantitative comparison of the approximation formula ( |3.8| ) with numerical results for 
acceptance rates justified our approximations in three different models. For the remainder of 
this section we assume that the approximation formula is valid. We now discuss the kinematical 
behavior of two other spin models that are asymptotically free in two dimensions. 



6.4 The two dimensional 0{N) nonlinear cr-model 

The two dimensional 0{N) nonlinear a-model is defined by the partition function 



Z = Jl[{d''s,Sisl-l))eMP E ^ 



(6.14) 



<x,y> 



where the Sx are iV-component real unit vector spin variables. The case N = 2 corresponds to 
the XY model. For N > 2 the 0{N) nonlinear a-model is in a disordered phase for all values of 
f3. In the limit (3 -^ oo, the correlation length diverges exponentially in f3. The critical behavior 
of the model is described the perturbative renormalization group which implies asymptotic 
freedom. 

Nonlocal Metropolis updates can be defined in close analogy to the XY model by updating 
in f/(l) subgroups of 0{N). Let us consider the special case where the updates act on the first 
two components of the s^ variables by 



Rfr 



(6.15) 



with 



e(3) 



v 



,{N) 



Rx 



sin(s'?/', 




V 







- sm{silJx) 

cos^sipx) 

1 





o\ 





/ 



(6.16) 



Again ip obeys the normalization condition ( p.8|) . Other directions of the U{1) subgroup in 
0{N) can be obtained by applying a randomly chosen global 0{N) rotation to the configuration 
before updating according to eq. ( |6.15| ). 
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The nonlocal updating proposal ( 6.15|) changes the Hamiltonian by 



A7^ = -P Yl {(^xSx) ■ {RySy) - Sx ■ Sy} 



<x,y> 



-^ E {4'^4'Mcos(.(^.-^,))-l] + .«3f sin(.(^.-^,)) 



<x,y> 



-.(^).«sin(. 



(^x - ^y)) + Sfhf^ [C0S(S(^. - ^,)) - 1]} . (6.17) 

The Hamiltonian of the 0{N) model is invariant under a global reflection of a fixed component 
sW —^ —s^\ Therefore all expectation values which change sign under such a transformation 
vanish. Using this, we obtain for hi = (A7i) 

hi = PE2 E [1 - cos(s(^, - ijy))] , (6.18) 

<x,y> 

with 

^2 = (4')4'^ + 4')4'^) = |(^.-^,), (6.19) 

where x and y are nearest neighbors. For piecewise constant kernels, hi is proportional to Lb- 
For smooth kernels hi will become independent of Lb for large enough blocks. For small s, 

hi ^ y^(3E2 J2 (^- - ^y? = b^(^E2 a , (6.20) 

<x,y> 

with a = {ip, —Alp). This quantity becomes independent of Lb for large block sizes. 

As in the XY model, the scale dependence of acceptance rates in the 0{N) nonlinear 
o"-model with A^ > 2 is like in massless free field theory. This behavior of the acceptance rates 
was already observed numerically in multigrid Monte Carlo simulations of the 0(3) nonlinear 
(T-model in two dimensions with smooth and piecewise constant kernels |T^ . Using a V-cycle and 



smooth interpolation in a unigrid implementation, a dynamical critical exponent z ~ 0.2 was 
observed. With a V-cycle and piecewise constant interpolation the result was z ^ 1.3. These 
values for z are very similar to the results for z in massless free field theory (cf. section ^). 



In the 0(4) model, z ^ 0.6 was found |jT9|, using piecewise constant interpolation and a 
W-cycle in a recursive multigrid scheme. 

One can do an analogous discussion for nonlinear cr-models with global SU{N) x SU{N) 
invariance, leading to a similar prediction for the scale dependence of the acceptance rates. In 
numerical experiments with a multigrid algorithm for the SU{3) x SU{3) nonlinear cr-model 



z ~ 0.2 was found 16 . 



6.5 The two dimensional CP^ ^ model 

A different class of nonlinear a-models in two dimensions are the CP'^^^ models. The complex 
projective space CP^~^ is a complex manifold with the topological structure of a 2A^ — 1 
dimensional sphere, where points related by a U{1) transformation are identified. A simple 
realization of a lattice action with a local U{1) invariance is 

n{z) = 2(3 Y. [^-\'^-^y?] (6.21) 

<x,y> 
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where the z^ are iV-component complex unit vectors. The CP"^ model is equivalent to the 0(3) 
nonlinear cr-model. For all A^, the CP^~^ models are in a disordered phase at all values of /3. 
As for all nonlinear cr-models in two dimensions, the critical behavior of the model is described 
by the perturbative renormalization group which implies asymptotic freedom for /3 — > cxd. In 
this limit the correlation length diverges exponentially in /3. 

Nonlocal Metropolis updates are formulated by updating in an embedded two dimensional 
XY model. Consider the special case where the updates act on the real and imaginary part of 
the first component of the z^. variables by 



^X^X 1 



with 



(6.22) 



/ 4^^ \ 

'^x 



Rx 



( exp(is^2,) 




V 









1 









1/ 



(6.23) 



As above, ijj is normalized according to eq. ( p.8|) . Other orientations of the embedded f/(l) sub- 
group can be obtained by acting on other components of the ^^.-variables, where the component 
is selected randomly. 

The energy change of such an update is 



AH = -2(3 Y. \^RxZx-R 



<x,y> 



.yZy\ 



N 



I Zx ' Zy I 



-4/3 E Re z«z«E4 

<x,y> I i=2 

■)N-1 



).W 



-is{i)x-i)y) _ j^ 



(6.24) 



The Hamiltonian of the CP model is invariant under a global complex conjugation z^ ^^ z^. 
Therefore the expectation value 



Q 



.(1)^(1) 



N 

E 

i=2 



(0 .^) 



Z^ Z- 



(6.25) 



where x and y are nearest neighbors, is real. Using this property, we find for the average energy 
change of the update ( |6.22| ) 



/ii = 4/?Q Y. [l-cos(s(^,.-V,))] . 

<a;,j/> 



(6.26) 



For piecewise constant kernels, h\ is proportional to L^- For smooth kernels hi will become 
independent of L^ for large enough blocks. For small s. 



hi ^ 2s^(3Q Y. (^- - '4^yf = 2s'/?g « • 

<x,y> 



(6.27) 



As above, a = (■?/', —Aip). a becomes independent of Lb for large block sizes. 

As in the XY model and in the 0{N) nonlinear cr-model, the scale dependence of acceptance 
rates in the CP^~^ model is like in massless free field theory. 

In multigrid Monte Carlo simulations of the CP^ model with smooth interpolation and a 
l^-cycle in a unigrid implementation, a dynamical critical exponent z ~ 0.2 was found [|T6[. 
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6.6 Summary 

A detailed quantitative analysis of the scale dependence of Metropolis updates was performed. 
Our approximation formula has proven to be quite precise in three different models. The results 
for all models are consistent with the following rule: 

SufEciently high acceptance rates for a complete ehmination of CSD can only be 
expected if hi = {7i{(f) + stp) — ?i(0)) contains no relevant operator in ip. 

As we have seen above, the typical candidate for a relevant operator in ip is always an 
algorithmic mass term of the type ~ s"^ J2x '4^1 with degree of relevance r = 2. 

This rule is formulated for smooth kernels. For piecewise constant kernels, it has to be 
modified. There, the kinetic term a oc Lb is relevant as well. At least in free field theory this 
disadvantage can be compensated for by using a W-cycle. Apart from this modification the 
rule carries over to the case of piecewise constant kernels. 

With the help of this rule it is possible to predict whether a given method will have a chance 
to eliminate CSD. 

Let us reformulate the rule in heuristic way. Recall the normalization condition ( p.8| ) for 
the interpolation kernels ip 

^B ^B Z_^ Wx = Ox',x'^ ■ 

Xdx' 

The factor L^ in front of the average L^ Hxex'i- ■ ■) was chosen such that the behavior 

of different terms in ip was analogous to the perturbative renormalization group. We observed 
that 

T.€ - Ll (6.28) 

X 

for piecewise constant and for smooth interpolation. The quantity a behaved like 

a = {ijj, -Aijj) ~ const (6.29) 

for smooth interpolation and 

a = {ip, -A^) ~ Lb (6.30) 

for piecewise constant interpolation, independent of the dimension d. See also the definition of 
the degree of relevance at the end of section ^. 

Let us switch to dimensionless kernels 

(6.31) 











i^x 


-Lf 


-2)/2 - 


which 


are 


now 


normalized according 


to 

T-d 

^B 


>;4 

x&x' 


= 5x',x'„ ■ 



(6.32) 
E.g. a d-dimensional piecewise constant interpolation kernel in this normalization is just 
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Then we find 

E^x ~ Li, (6.34) 

X 

for piecewise constant and for smooth interpolation. A kinetic term in the kernels ip ^low 
behaves as 

a = {ip, -AiP) ~ L|-2 (6.35) 

for smooth interpolation and 

a = (^, -AiP) ~ Lf^ (6.36) 

for piecewise constant interpolation. 

Assume that we use piecewise constant interpolation kernels ip. Then the terms that appear 
in the average energy change hi = {H{4> + sip) — 7i(0)) of a nonlocal update can be interpreted 
as follows: 

• a kinetic term {ip, —Aip) is a surface term, i. e. the energy costs of a nonlocal update are 
proportional to the surface of the block 

• a mass term ~ J2x W is a volume term,, i. e. the energy costs of a nonlocal update are 
proportional to the volume of the block 

There are two different ways how multigrid algorithms can deal with a surface term: First, one 
can put more emphasis on the updating of coarse levels by using a W-cycle. Second, one can 
switch from piecewise constant interpolation to smooth interpolation. Then the energy change 
at the surface of a block is spread uniformly over the entire block, leading to {ip, —Aip) ~ Lg . 

However, up to now there is no way how to overcome the effect of a volume term that causes 
CSD in multigrid Monte Carlo simulations of critical models. 

The different role of surface and volume terms seems to be a general feature of nonlocal 
update methods that multigrid algorithms share with cluster algorithms. In cluster algorithms 
nonlocal moves are performed by stochastically growing a cluster of spins and then flipping all 
spins in the cluster simultaneously. In ref. it is argued that if the energy cost of this global 
flip is a "bulk term" that is proportional to the volume of the cluster, a cluster algorithm is 
expected to work badly. 

The natural way to construct nonlocal piecewise constant updates that only have surface en- 
ergy costs is to look for a global symmetry of the Hamiltonian of the model under consideration. 
Then the idea is to apply the corresponding symmetry transformation locally on blocks. E. g. 
in the massless free field theory, the global symmetry is the shift symmetry (px —>■ (px + const, 
which is then applied locally on blocks. In the XF-model, the global (9(2)-symmetry of the 
model is applied locally by performing piecewise constant 0(2)-rotations on blocks. 

To summarize this discussion we give an intuitive Leitmotiv for the development of new 
multigrid methods: 

A piecewise constant update of a nonlocal domain should have energy costs propor- 
tional to the surface of the domain, but not energy costs proportional to the volume 
of the domain. 

We will use this heuristic criterion in the development of multigrid algorithms for nonabelian 
gauge theories. The intuitive Leitmotiv can then be made precise by performing a quantitative 
acceptance analysis. 
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7 Multigrid Simulation of the Sine Gordon model 

7.1 Motivation 

In this section we report on a multigrid Monte Carlo study of the Sine Gordon model in two 
dimensions. We have two questions in mind: 

First, we want to check the theoretical prediction that a W-cycle (cycle control parameter 
7 = 2) with piecewise constant interpolation will not eliminate CSD in the rough (massless) 
phase of the Sine Gordon model. Recall the analysis of this model in the previous section: The 
Sine Gordon model was an example for a critical model where the expansion of (7i(0 + ip)) 
contained a relevant (mass) term. If such a term is present, Metropolis step sizes e{LB) on 
block lattices with increasing block size Lb have to be scaled down like e{LB) ~ 1/-^b in order 
to obtain block size independent acceptance rates. This strong decrease of step sizes on large 
blocks will lead to CSD. 

Second, we want to ask the question whether one can circumvent this slowing down caused 
by too small steps on large blocks by accumulating many of these steps randomly. 

We applied an intuitive random walk argument for free field theory in section ^. There 
this reasoning explained why one can compensate for the disadvantage of £:(Lb) ~ I/^/Lb with 
piecewise constant interpolation by using a W-cycle with 7 = 2. We can use the same argument 
when the Metropolis step sizes decrease even stronger: 

A constant accumulated step size on all length scales could in principle be achieved in the 
following way: For step sizes scahng down like s^Lb) ~ ^/ Lb and a coarsening by a factor 
of two, the Metropolis step size on a next coarser grid is too small by a factor of two. If we 
assume that subsequent update steps within a multigrid cycle are statistically independent and 
that these steps accumulate in a random-walk like way, we can expect to compensate for this 
decrease by increasing the number of updates on the next coarser grid by a factor of four. This 
can be achieved by a higher cycle with cycle control parameter 7 = 4. By applying a higher 




Lb = 1 
Lb = 2 
Lb =4: 
Lb = S 



Figure 8: Higher cycle with 7 = 4 

cycle with 7 > 1, 7 times more updates on each coarser lattice are performed. An example for 
a higher cycle with 7 = 4 is given in figure §. 

For a recursive multigrid algorithm, the computational effort is ~ L'^ for 7 < 2^^ and 
~ L'^ log L for 7 = 2*^ in (i dimensions j^lj . Therefore a higher cycle with 7 = 4 is practical for 
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d > 2 and borderline practical ioi d = 2. 

In summary, if the random- walk picture is correct, we can expect to overcome the problem 
of strongly decreasing step sizes on large scales by doing more work on coarser lattices. 

7.2 The multigrid algorithm 

The Hamiltonian of the two dimensional Sine Gordon model is 

n{ct>) = ^(0,-A0) -CEcos27r0, (7.1) 

on an L X L lattice Aq. Our simulations are organized as follows: In order to allow for high cycle 
control parameters 7, we use a recursive multigrid algorithm, piecewise constant interpolation 
and a coarsening with a factor of two, as illustrated in figure [^. Let us describe the details of 
the recursive multigrid algorithm. An update to the current fine grid field (px is given by 



$^/ iixex' . (7.2) 



The point x on the fine grid is contained in the 2x2 block x' on the coarse grid, and $3,/ 
denotes the displacement fieldQ corresponding to the block x'. 

The conditional coarse grid Hamiltonian is calculated by substituting eq. jfT^ ) into eq. ( [7.1|) . 
Its general form on the block lattice A^ is 

^fc(*) = f Y. i'^x'-%'f- Y. fx'^^x'-C Y (ax'Cos27r$^. + 6^.sin27r$^/) . (7.3) 

(x',j/'>GAfc x'eAfe x'eAfc 

Note that the fundamental Hamiltonian Ti on the fine grid Ag is a special case of this general 
coarse grid Hamiltonian Ti^ on A^ with 

(7.4) 
(7.5) 
(7.6) 
(7.7) 

The set of coefficients {a', /', a', b'} of TCk+i on the next coarser grid Afc+i can be calculated 
from the set of coefficients {a, /, a, b} on A^t by recursion: 

a' = 2a , (7.8) 



a = 


= i/P, 


f. - 


= for all X E Aq , 


O'x - 


= 1 for all X G Ao , 


K -- 


= for all X e Aq . 



fL' = E 



f.-2aJ2 



■'X 



x&x' L yn.n.x 



(7.9) 



«L' = E (^^^^os^vr^^. + 6x'Sin27r0^.) , (7.10) 

K' = E (~"^si^27r0^ + &xcos27r</)^.) . (7.11) 

x^x' 

Here (px denotes the configuration on A^ when we switch to the next coarser grid. 



^Although we use the same letter $, the displacement field as used here and the block spin field as introduced 
in section differ by a constant 
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For the updating of the variables on A^ within the multigrid cycle, we choose a sweep of 
single hit Metropolis updates before the coarsening and after the interpolation, respectively. 
The maximum Metropolis step size s{Lb) is scaled down like 1/Lb- Then, acceptance rates 
of approximately 50% are observed on all block lattices, in accordance with the theoretical 
analysis of section ^. We tested the implementation of the algorithm by comparison with exact 
results in the Gaussian model {( = 0) and against Monte Carlo data obtained on small lattices 
by a cluster algorithm for the Sine Gordon model [H^ . 

7.3 Simulation and results 

We study the dynamical critical behavior of the different versions of the algorithm in the rough 
phase, where the correlation length is infinite and the physical length scale is set by the linear 
size of the lattice L. Thus, we expect the autocorrelation time r to diverge with the dynamical 
critical exponent z like r r^ L^. 

As observables we take the energy 

E = ^^T.i<t>.-<^yf (7.12) 



and the surface thickness a^ as defined in eq. (|6.2|). 

The simulations were performed aX j3 = 1.0, C, = 0.5. This is deep in the rough phase. Note 
that in the limit C ^ oo which corresponds to the discrete Gaussian model the critical coupling 
is j3c = 0.7524(8) [^. Starting from an ordered configuration, measurements were taken after 



equilibration at each visit of the finest lattice. The run parameters and results are summarized 
in table ^ for the W-cycle and in table ^ for the higher cycle with 7 = 4. 

From the autocorrelation function for the observable A 

/,N _ {AsAs+t) - {A) . . 

PAit) - ^^2) _ ^^)2 ' (7-13) 

we compute the integrated autocorrelation time 

-| oo 

r^nt,A = ^ + EPaW- (7-14) 

^ t=l 

The asymptotic decay of the autocorrelation function pA{t) is described by the exponential 
autocorrelation time Texp,A'- 

PA(t) ~ exp{-t/Texp,A) for t ^ oo . (7.15) 

If pA{t) is a single exponential, Tint,A ~ ^exp.A- 

The errors on Tint,A were calculated by a window method |48|. For the self consistent 
truncation window we took 4x^4, which is sufficiently large if the autocorrelation function 
shows an exponential decay. This is indeed the case, see figure |[ 

The numerical results for the autocorrelation times Tint are listed in table ^ for 7 = 2 and in 
table ^ for 7 = 4. The values for Tint are measured in the number of visits on the finest lattice. 
Note that our runs are longer than 10 000 Tint (longer than 4 000 Tint on the 128^ lattice). Figure 
pni shows the dependence of Ti^t^a'^ on L for 7 = 2 and 7 = 4. For comparison, we plotted lines 
which correspond io z = 2. 
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Table 3: Numerical results for the W-cycle (7 = 2) in the 2-d Sine Gordon model on L x L 
lattices in the rough phase, P = 1.0, ( = 0.5 



L 


statistics 


discarded 


E 


Tint,E 


a' 


'^int,a'^ 


4 


25 000 


2 000 


0.934(4) 


0.90(3) 


0.268(1) 


0.96(3) 


8 


50 000 


2 000 


0.986(1) 


0.97(2) 


0.3809(9) 


1.35(3) 


16 


100 000 


2 000 


0.9956(5) 


1.04(2) 


0.4896(7) 


2.70(8) 


32 


300 000 


2 000 


0.9987(2) 


1.03(1) 


0.5996(7) 


8.54(19) 


64 


500 000 


2 000 


0.99945(5) 


1.04(1) 


0.7105(10) 


30.5(1.0) 


128 


500 000 


4000 


0.99966(3) 


1.04(1) 


0.8218(19) 


113.7(6.9) 



Table 4: Numerical results for the higher cycle with 7 = 4 in the 2-d Sine Gordon model on 
L X L lattices in the rough phase, j3 = 1.0, ( = 0.5 



L 


statistics 


discarded 


E 


Tint,E 


a' 


'^int,a'^ 


4 


25 000 


2 000 


0.940(4) 


0.89(3) 


0.268(1) 


0.91(3) 


8 


25 000 


2 000 


0.986(2) 


0.94(3) 


0.380(1) 


1.14(4) 


16 


25 000 


2 000 


0.9965(10) 


0.94(3) 


0.488(1) 


1.67(6) 


32 


100 000 


2 000 


0.9985(3) 


0.95(1) 


0.5997(9) 


4.15(11) 


64 


300 000 


2 000 


0.99945(7) 


0.96(1) 


0.7113(9) 


14.2(4) 


128 


300 000 


2 000 


0.99962(4) 


0.95(1) 


0.8213(18) 


58.2(3.3) 
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If we fit our data for the autocorrelation time of the surface thickness in the range 32 < L < 
128 with the Ansatz Tint,a^ = cL^, we obtain z = 1.86(4) with x^/dof = 0.22 for the W-cycle 
(7 = 2), and z = 1.86(4) with x^/dof = 4.6 for the higher cycle (7 = 4). The uncertainty in 
z is dominated by the relative error of Tint on the largest lattice, which is about 6% in both 
cases. We therefore estimate 

z = 1.9(1) for the W-cycle (7 = 2), 
z = 1.9(1) for the higher cycle (7 = 4). 



7.4 Summary 

In this section we reported on a simulation of the two dimensional Sine Gordon model with a 
multigrid Monte Carlo algorithm. Using a recursive multigrid implementation with piecewise 
constant interpolation we studied the dynamical critical behavior of the two variants of the 
algorithm (7 = 2 and 7 = 4) in the rough phase. From this study we confirm that, as 
already predicted in section ^, CSD is not reduced by a W-cycle with piecewise constant 
interpolation. According to the acceptance analysis, we would expect the same result with 
smooth interpolation. Moreover, the results clearly show that the attempt to compensate for 
decreasing acceptance rates on large blocks by choosing a higher cycle with 7 = 4 does not 
improve the dynamical critical behavior of the algorithm. We conclude that a random-walk 
argument as stated above is not correct in the case of the Sine Gordon model. 
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Figure 9: Autocorrelation functions p{t) for the surface thickness a"^ in the rough 
phase of the two dimensional Sine Gordon model, j3 = 1.0, ( = 0.5. Top: 64^ 
lattice, bottom: 128^ lattice. p{t) is shown for the W-cycle (7 = 2) and the higher 
cycle ('J = 4) respectively. 
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Figure 10: Dependence of the integrated autocorrelation time for the surface 
thickness a"^ on the lattice size L in the rough phase of the 2-d Sine Gordon 
model, j3 = 1.0, ( = 0.5. Errors are smaller than the symbols used. The lines 
correspond to z = 2. 
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8 Multigrid algorithms for lattice gauge fields in two di- 
mensions 

In this section, we start the discussion of multigrid algorithms for lattice gauge theory. We 
begin with a review of the method by Laursen, Smit and Vink for abelian gauge fields in two 
dimensions. An new difficulty is to change from an abelian to a nonabelian gauge group. The 
concepts that are needed for nonlocal updates of nonabelian gauge fields in in two dimensional 
SU{2) lattice gauge theory are discussed. We introduce a new updating procedure, the time 
slice blocking, and analyze its kinematics. 

Let us introduce the notations for lattice gauge theory in d dimensions. We consider parti- 
tion functions 

Z= f YldU^,, expi-mU)) . (8.1) 

The link variables [4.^^ take values in the gauge group f/(l) or SU{N), and dU denotes the 
corresponding invariant Haar measure. The standard Wilson action 'H{U) is given by 

'HiU) = /3E[1- ^ReTrf/p] . (8.2) 

V 

The sum in (|8.2|) is over all plaquettes in the lattice. The U-p are path ordered products around 
plaquettes "P, 

U* denotes the hermitean conjugate (= inverse) of U. 

8.1 Abelian gauge fields: U{1) in two dimensions 

8.1.1 The algorithm of Laursen, Smit and Vink 

A first multigrid Monte Carlo algorithm for abelian gauge fields was implemented and tested by 
Laursen, Smit and Vink for f/(l) lattice gauge theory in two dimensions [|18|. It was inspired by 
a deterministic multigrid procedure for the minimization of a Hamiltonian of the form arising 
in lattice gauge theory [^ . As an example and a starting point for the following discussion of 
nonabelian gauge fields we review their algorithm in the unigrid language. 

In the case of gauge group U{1) the link variables are parameterized with an angle — tt < 
6j:,f, < IT through 

Ux,f, = exp(2e^,^) . (8.4) 

Nonlocal updates are defined as illustrated in figure |TT]: One chooses a square block x'^ 
of size L^ and a direction r with r = 1 or 2. During the update, r will be kept fixed. All 
the link variables Ux^t attached to sites x inside the block x'^ are proposed to be "rotated" 
simultaneously: 

Ux,r -^ exp{istlJx)Ux,r , (8.5) 

or in terms of the link angles 

Ox,r ^ Ox,r + SlIJx • (8.6) 
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Figure 11: Blocking of link variables for two dimensional U{1) gauge fields 



ip denotes an interpolation kernel as introduced in section 
dimensional piecewise constant kernel 



t 



const 



1 for X G Xq 
for X ^ Xq 



In the reported work a two 



(8.7) 



was taken. For this choice of interpolation, a recursive multigrid implementation with a W- 
cycle is feasible [|l8l . Another interesting feature of this implementation is (in contrast to the 
time slice blocking below) that both directions of link variables pointing from the block x^ in 
the r = 1 and in the r = 2 direction can be updated on a coarse layer without going back to 
the finest lattice. 

The reported result for the dynamical critical exponent z for the simulation with a W-cjcle 
is z ~ 0.2 with respect to the autocorrelation times of Polyakov loop observables. However, it 
was also observed that the multigrid algorithm is not able to move efficiently between different 
topological sectors. This leads to autocorrelation times r ~ exp(c/3) with c ~ 1.9 for the 
topological charge. If the physical size of the lattice is fixed, this will lead to r ~ exp(c'L), 
i.e. super-CSD increasing exponentially in L. The exponent quoted above should therefore be 
interpreted with some care. 

Let us analyze the kinematical behavior of this algorithm with the help of eq. (|3.8|) 



n{s) 



erfc ( ^ 



For hi = (A7i) we find in the case of piecewise constant interpolation 

hi = 2(3PLb[1-cos{s)] = (3PLbs'^ + 0{s^) 



(8.8) 



with P = (Trf/p). The derivation of this expression is described in appendix ^. In order to 
obtain acceptance rates independent of the block size, we have to rescale the amplitudes for 
the nonlocal updates like s ~ 1/\/Lb, as in massless free field theory with piecewise constant 
interpolation. With respect to this, the behavior of this multigrid algorithm in the U{1) lattice 
gauge theory in two dimensions is similar to a multigrid algorithm in massless free field theory. 
Therefore the reported acceleration for Polyakov loop observables can be understood by our 
analysis. 
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Figure 12: Decoupling of time slices A[ and Aj . 

On the other hand the problems with topological quantities are again an example (cf. the 
discussion of the two dimensional XF-model in the vortex phase in section |]) for the fact that 
good acceptance rates alone are not sufficient to overcome CSD. 



8.1.2 Comments and possible modifications 

Let us denote the time slice of lattice sites with r-component t as A[ = {x G Aq | Xt- = t}. 
Here the name "time" direction has no physical meaning. We use this word to label the fixed 
direction of link variables that are updated simultaneously. Of course, the time direction in an 
update algorithm will be changed periodically from r = 1 to r = 2. In the following, we will 
denote the time direction with r and the spatial direction(s) different from r with /i. 

In the unigrid picture a general feature of nonlocal updating schemes in lattice gauge theory 
is transparent: As long as we restrict the possible nonlocal changes in the configuration to link 
variables Ux^t of a fixed time direction r and keep all other variables f/^^^ with fi y^ t unchanged. 



adjacent time slices decouple. This is illustrated in figure ^. 

Here, link variables Ux^t pointing from sites in two different adjacent time slices A[ and A2 
are shown. The point is that there are only plaquette terms in the Hamiltonian that contain 
link variables Ux,t pointing from sites x in the same time slice A[, i.e. either in A[ or A^. This 
means that the link variables pointing from sites in A[ and from sites in Aj are independent as 
long as all other U^^fj, with fi ^ t are fixed. 

A possible modification of the algorithm is the use of smooth interpolation kernels ip. If we 
repeat the acceptance analysis with the help of formula (|3.8| ) for a general kernel ip we obtain 
for hi = (A7i) (cf. appendix 0) 



hi = PP Yl [1 - COs(s(^a;+A " V'x))] 



^.9) 



a; G Ao 



The sum does not include contributions from links which point in the time-direction r. There- 
fore hi is a sum of independent contributions from time slices AJ" orthogonal to the r-direction. 
This is a consequence of the independence of time slices. In particular, no smoothness of kernels 
is needed in the r-direction. From now on we choose ipx to be constant in that direction and 
assume that the support of ip in r-direction is restricted to the Lb time slices A[ that contain 
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Figure 13: Geometry of the time slice blocking: The marked link variables are updated 
simultaneously. The bottom of the block is indicated by a dashed line. 
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(8.11) 



with «! = {ip', —Alp'). Here, ip' denotes the kernel ip restricted to the one dimensional time slice 
AJ". A factor of L^ is absorbed in ip' such that ip' is properly normalized as a one dimensional 
kernel. 

This means that for smooth interpolation we can obtain acceptance rates independent of the 
block size by choosing the amplitudes for the nonlocal updates like s ~ const. The kinematical 
behavior of this version of the multigrid algorithm is similar to a multigrid algorithm with 
smooth interpolation in massless free field theory. Therefore we expect a successful acceleration 
for Polyakov loop observables also for this version. 

In summary, the independence property of time slices has several consequences: First, 
interpolation kernels ip do not need to be smooth in the time direction. Second, since updates 
of link variables in different time slices are statistically independent, two different nonlocal 
update schemes are possible: the time slice blocking and the square blocking. In the time slice 
blocking method updates of one dimensional blocks of size Lb are performed on separate time 
slices A[ in sequence. In the square blocking scheme as used by Laursen, Smit and Vink one 
builds Lb x Lb blocks out of "staples" of Lb one dimensional blocks of size Lb and performs 
the updates on this square block simultaneously. The analysis of the kinematics is the same 
for both schemes. For simplicity we are going to adopt the time slice blocking in the following. 

An important point is that the decoupling of time slices is independent of the gauge group 
and carries over to higher dimensions. We are going to use this fact as a basic ingredient of 
nonlocal updating schemes for nonabelian gauge fields. 

8.2 Nonabelian gauge fields: SU{2) in two dimensions 
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Figure 14: Change of plaquette during naive nonlocal updating 



8.2.1 The nonabelian character of the gauge field 

Our method for a nonlocal updating procedure for nonabelian SU{2) gauge fields in two dimen- 
sions is based on the time slice blocking. It is illustrated in figure |T^. We start the discussion 
with a naive generalization of the updates in the abelian case: Choose a one dimensional block 
x'^ of size Lb within a time slice A[ and update all link variables Ux^r pointing from this block 
in the r-direction 



U.. 



Ul 



R'xUxT for all X E x' 



where the "rotation" matrices Rx are in SU{2). We parametrize them as 

Rx{n, s) = cos(s^a;/2) + isin^sipx/^) n-a , 



(8.12) 



(8.13) 



where n denotes a three dimensional real unit vector, and the at are Pauli matrices, n will be 
taken randomly from the three dimensional unit sphere, and ip will have support on the one 
dimensional block x'^. 

The simplest version is a piecewise constant "rotation" where Rx is a rotation matrix R 
independent of x. Let us examine how a plaquette in the interior of the block (as illustrated in 
figure |l^ changes under this update: 



Fold = Yli{UiU2U-*Ul 
1 



\TY{UrRU2U^{RUi 



^ new 2 — \- J^ ^-c\ •*/ / r) 

with the notations as in figure 0. Let denote the links 

(x, X + ii) for all X, x + /i G x'^ 



Tt{R*UiRU2U;U^) 



(8.14) 



(8.15) 



as the links in the bottom of the block, as marked in figure |T3|. In the abelian case the rotation 
matrix R and the link variable in the bottom of the block Ui would commute, and we would 
have Fold = Pnew, i-e. the plaquettes in the interior of the block would remain unchanged. 

But since our gauge group is nonabelian the resulting change of the plaquette by this naive 
generalization would lead to energy costs proportional to the volume of the block. 

Recall the intuitive Leitmotiv for the development of new multigrid methods from section |^: 
A piecewise constant update of a nonlocal domain should have energy costs proportional to the 
surface of the domain, but not energy costs proportional to the volume of the domain. 

Energy costs proportional to the volume of the block lead to an algorithmic mass term 
that suppresses the amplitudes of the nonlocal moves on large scales. Therefore the naive 
generalization will suffer from CSD. 
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This does not come as a surprise. Due to the gauge invariance of the model, the hnk variables 
do not have a gauge invariant meaning. Therefore the rotation R that is constant over the block 
for a given gauge can be arbitrarily rough and disordered after a gauge transformation. It is 
therefore natural to assume that the rotation matrices have to be chosen in a gauge covariant 
way. 

8.2.2 Gauge covariant time slice blocking algorithm 



We use the additional gauge degrees of freedom and generalize (|8.12) to 



U,,r -^ K,r = RMU.,r for all X G < , (8.16) 

with Rx{g) = QxRxQx and (/-matrices g^ G SU{2). In the abelian case we obtain nothing new, 
because g^ and Rx commute. In the nonabelian case a plaquette in the interior of the block 
changes under a piecewise constant rotation R according to 

Fold = Yli{UrU2U;Ul) -. P„e«, = iTr(i?*[/fi?[/|f/|*f/f ) . (8.17) 

The gauge transformed link variables are given by U^^^ = gxUx,^j.gx+^- If we choose the 
5f-matrices in such a way that the link variable Uf in the bottom of the block is equal to 
unity, it will commute with an arbitrary rotation matrix R. Then the plaquettes in the interior 
of the block stay unchanged as in the abelian case. 

From this discussion we are led to the following gauge condition: Choose the (/-matrices in 
the bottom of the one dimensional block x'^ such that 

Ul^ = 9xUx,^.gUf, = l for all {x,x + fi) e x'^ . (8.18) 

(The bottom of the block is indicated in figure |13|.) In other words: All gauge transformed link 
variables in the one dimensional bottom of the block should be equal to unity. We denote this 
gauge condition as block axial gauge. Note that we still have the freedom of a constant gauge 
transformation within each block x', 

gx —>■ hx'gx for all a; G x' . (8.19) 

Although we use the term "gauge condition" here, we do not intend to perform the gauge 
transformation g. We use the concept of gauging only to define covariant rotations Rx{g) = 

glRxgx- 

The gauge transformation properties of these updates are as follows: If we apply an arbitrary 
gauge transformation h to the gauge field U 

Ux,, - Ul^ = hxUx,,h:^i, (8.20) 

the (/-matrices transform like 

9x -^ g^x = 9xK (8-21) 

and the covariant rotation matrix Rx{g) = g*xRxgx transforms according to the adjoint repre- 
sentation: 

Rx{9) ^ {Rx{g)f = hxRx{g)K ■ (8-22) 

As a consequence, if we apply the updates to the gauge transformed configuration U^ 

f/,V - Kr)' = iRMfU^x,r = hxRx{g)KhxUx,rhUr 
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= KRMU.,rK_,, = KK^^K^, = (f/;,)'^ , (8.23) 

the updating commutes with a gauge transformation h and is therefore gauge covariant. 

Let us summarize the steps of the time shce blocking scheme for SU{2) in the unigrid 
language: 

1. Choose a block x'^ of size Lb that is contained in the slice A[. All link variables Ux^t 
pointing from sites x inside the block in r-direction will be moved simultaneously. 

2. Find the gauge transformation g defined by the block axial gauge condition ( |8.18| ) such 
that f/^ = 1 for all link variables in the bottom of the block. 

3. Propose new link variables U'^ ^ by 

Ux,r^K,r=RA9)Ux,r, (8.24) 

with Rx{g) = glRxQx and 

Rx{n, s) = cos(s-?/'s/2) + i sm{sipx/'2) n-a . (8.25) 

s is a uniformly distributed random number from the interval [— e, e], n is a vector selected 
randomly from the three dimensional unit sphere, and ?/; is a one dimensional kernel. 

4. Calculate the associated change of the Hamiltonian A7i and accept the proposed link 
variables with probability min[l,exp(— A7i)]. 

The detailed balance condition is fulfilled by this updating scheme: For the naive version 
with (7 = 1 it is straightforward to show that the detailed balance condition holds, since the 
rotation matrices Rx are chosen according to a probability distribution which is symmetric 
around unity. 

If we now take g according to some gauge condition, we have to be careful that we get the 
same g before and after the move Ux^t ~^ f^x t i^ performed. Otherwise this move would not 
be reversible. In other words: The gauge condition yielding g must not depend on Ux,t- This 
is indeed the case, since only link variables Ux,^i with \i ^ r enter in the block axial gauge 
condition. 

The details of an implementation and simulation of the covariant time slice blocking algo- 
rithm for SU{2) gauge fields will be described in section ^ 

8.2.3 Acceptance analysis of the proposal 

The energy change associated with the update proposal ( p.l6| ) is 



A^ = -f E Tr(f/^ - Uv) = -f E Tr{{Rx{grUx,,Rig)x+f. - Ux,,)H:J , (8.26) 

where H*^^ = Ux+ji^tU*^^ fj_U*^^, and ip stands for a one dimensional interpolation kernel. The 

relevant quantity for the acceptance rates is hi = (A7Y). If we assume that the (yf-matrices are 

chosen according to the block axial gauge condition and that ip vanishes outside the block x'^ 

we find 

hi = (3P E [1-cos(s(7/>,-V'x+a)/2)] 
xeAi 
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2 2 

= '-pP Y. {^.-i^.+^f + 0{s') = '-pPa, + 0{s') (8.27) 

with P = (iTrt/p) and ai = (^z^, — A^). The derivation of eq. (|8.27|) is described in appendix 0. 



Remember that for the time shce blocking in two dimensions ijj is a. one dimensional kernel. 

Thus the kinematical behavior of this method is the same as that of the massless Gaussian 
model in one dimension. In contrast to U{1) lattice gauge theory in two dimensions we do 
not expect any additional topological problems to occur for SU{2). Therefore we expect a 
successful acceleration of the simulation by the proposed time slice blocking algorithm. This 
prediction will be verified in in section |^. 

8.3 Summary 

In this section we discussed multigrid algorithms for lattice gauge theory in two dimensions. 
The successful acceleration of Polyakov loop observables in the simulation of abelian gauge 
fields in two dimensions by Laursen, Sniit and Vink can be understood by the kinematical 
analysis. 

An important observation is the statistical decoupling of adjacent time slices as long as only 
link variables in the time direction are updated. The time slice blocking algorithm is based on 
this property. The statistical independence of adjacent time slices is independent of the gauge 
group and of the dimensionality. 

The nonabehan character of the gauge field was discussed. We proposed the gauge covariant 
time slice blocking for SU{2) in two dimensions that is particularly adapted to the nonabelian 
character. From the kinematical analysis we expect a strong reduction of CSD in a simulation 
of the proposed algorithm. 
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9 Multigrid Monte Carlo simulation oi SU{2) lattice gauge 
fields in two dimensions 

In this section we simulate SU{2) lattice gauge theory in two dimensions by a multigrid algo- 
rithm. The acceptance analysis of the previous section predicted that CSD can be eliminated 
by the time slice blocking. To verify this statement, we perform numerical experiments on sys- 
tems with lattice sizes up to 256^. The details of the implementation of the time slice blocking 
are discussed. Then the run parameters and numerical results are reported. 

9.1 Implementation of the time slice blocking 

For a concrete implementation of the time slice blocking method as introduced in section |] we 
choose a unigrid algorithm with piecewise linear interpolation and a V-cycle. 

9.1.1 Smooth nonlocal heat bath updates in f/(l) subgroups 

Although the analysis of nonlocal updates in section ^ was performed in terms of a Metropolis 
version, we are going to use a heat bath version of the nonlocal moves in our simulation. We 
think that the change of the update method from Metropolis to heat bath will not affect the 
dynamical critical behavior in a substantial way. The advantage of the heat bath updating is 
that there are no tuneable parameters such as the Metropolis step size e{LB)- 

A heat bath implementation of nonlocal time slice blocking updates is possible if we use 
one dimensional piecewise linear interpolation and update in f/(l) subgroups of SU{2). This 
will be described in the following. 

Compared to section |^, we formulate the piecewise linear interpolation in a different lan- 
guage: Assume that the (yf-matrices defined according to the block axial gauge condition (|8.18|) 



have been applied as a gauge transformation in the bottom of the block x'^. Then the gauged 
link variables in the bottom of the block are equal to unity. Now a piecewise linear block update 
is formulated by multiplying the Lb gauged link variables ?7|^ pointing in r-direction from left 
to right by the S'f/(2)-matrices R, R^, i^^ . . . i?^^/^ R^b/2^ R^bI2-i^ . . . , i?. i? is given by 

R{n,e) = cos{e) + ism{e)n-a, (9.1) 

where the randomly chosen three-dimensional vector n specifies the direction of a U{1) subgroup 
in SU{2). All changes of plaquettes that are generated by this update can be written in the 
form 

Tr(f/^) = Ti{RV) = Tt{V) cos(^) + TT{in-aV) sin(^) , (9.2) 

with the SU{2) matrix V = Uj, or V = f/p*. By summing over all changed plaquettes ( pl2|) we 
obtain an overall change in the Hamiltonian of the form 

n{U') = a cos(e) + b sin(^) + const , (9.3) 

with real constants a and b. To generate U{1) random numbers distributed according to the 
distribution 

dprob(^) oc e'^'^os{0)+bsin{e)^Q ^g 4^, 

we use the fast vectorizable method of Hattori and Nakajima 

46 



9.1.2 Sequence of updates 

The sequence of updates is organized as follows: We start with a V-cycle of time slice blocking 
updates on the links Ux,t pointing in the r = 1 direction. The largest block size is L/2 on a L x L 
lattice. This means that the sequence of updated block sizes is Lb = 2,4,... L/2, L/2, L/4 ... 2. 
When all time slices have been updated by a V-cycle, we perform a sweep of local SU{2) heat 
bath updates through all links on the lattice. Here we use the "incomplete" Kennedy- Pendleton 
50| heat bath algorithm with one trial per link. This means that in a sweep through the lattice 



not all links but only a very high percentage of them are updated. The advantage is that scalar 



operations on a vector computer are avoided |51]. Then we do a V-cycle on all links pointing in 



the r = 2 direction and again a local heat bath sweep. This sequence is repeated periodically. 

Measurements are performed after each local heat bath sweep. To avoid effects from fixed 
block boundaries, we use stochastically overlapping blocks [|16] by applying a random translation 
before each V-cycle. 

9.1.3 Axial gauge 

In order to save computer time we use a slight modification of the gauge condition. Recall the 
block axial gauge ( fj.l(j| ): Take the (yf-matrices in the bottom of a block x'^ such that 

U^f, = 9xUx,,,9l+f, = 1 for all (x, x + fl) e x'^ . (9.5) 

We modify it to the axial gauge: Take the ^f-matrices in the bottom of a time slice A[ such 
that 

U^f, = gxU:,^^,gl^p^ = 1 for {x,x + jl)eKl, x^ = l,...L-l . (9.6) 

In other words: all spatial links but the last link are gauged to one within the bottom of a time 
slice. 

The computational advantage of this slight change in the gauge condition is that for the 
block axial gauge the gf-matrices have to be calculated for each block lattice individually. For 
the axial gauge the (yf-matrices are the same for all block lattices and they do not change during 
the updating on different block lattices with different Lb- Therefore we have to calculate them 
only once before performing the entire V-cycle. The gauge covariance properties of the update 
is not affected by this modification. 

9.2 Simulation and results 
9.2.1 Observables 

The observables measured are square Wilson loops 

iy(J,J) = (iTr(t/(C,,,)), (9.7) 

where U{Cjj) is the parallel transporter around a rectangular Wilson loop Cjj of size I x L 
On an L X L-lattice we measure 1^(1, 1), W{2, 2), W{4:, 4), W{8, 8), . . . , W {L/2', L/2). Another 
important class of quantities is built up from timelike Polyakov loops. A Polyakov loop at the 
one dimensional spatial point r is defined by 

L 

P, = lTrnt/(.,i),2. (9.8) 
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We measure the lattice averaged Polyakov loop 




and the lattice averaged Polyakov loop squared 



(9.9) 



P' 



■EPr 

r=l 



(9.10) 



9.2.2 Run parameters and results 

In order to investigate the dynamical critical behavior of the multigrid algorithm, we simulate 
a sequence of lattices with fixed physical size L ?» 10^, where the correlation length ^ is related 
to the string tension hj k hj ^ = Xj \fK. Then we have to use lattice sizes and /3 values such 
that L^//5 is constant. If we choose this large ratio of L/^, finite size effects are negligible. 
The detailed run parameters are given in table ^. The quoted correlation length is calculated 



by the exact solution [Q in the infinite volume limit. The main results of this solution are 
summarized in appendix |B|. We started our runs from ordered configurations with all link 
variables set equal to unity. After equilibration, measurements were taken after each local heat 
bath sweep through the lattice. 



Table 5: Run parameters for the multigrid Monte Carlo simulation of two dimensional SU{2) 
lattice gauge theory. 



p 


4 


16 


64 


256 


1024 


L 


16 


32 


64 


128 


256 


e 


1.55 


3.22 


6.51 


13.05 


26.12 



The static results of the simulation are given in table p. All our results for the Wilson loops 
are consistent with the exact solution ( [B.4| ) in the infinite volume limit. Here only results for 
iy(^, ^) and W^(f , f ) are quoted, which are the loops of about the size of a correlation length 
squared. We observed that these loop sizes have the largest autocorrelation times among the 
Wilson loops. 

In general we found a very fast decorrelation of subsequent configurations in the Markov 
chain. Typically the autocorrelation function p{t) dropped to zero within errors after 3 — 5 
measurements (see figure 0). Therefore it was impossible to look for an exponential regime 
in the decay of p(t). We tried to estimate integrated autocorrelation times Tmt with a self 
consistent truncation window of 4rj„t according to the method of ref. |48]. The results for the 
integrated autocorrelation times are given in table |^. Estimates for Tint are only quoted if we 
observed that p{t) was positive in the entire interval from t = to t = 4rj„t. Note that Tint is 
defined such that Tint = 0.5 in the case of complete decorrelation. All our runs are longer than 
30 000ri„t. 

In summary, all r's are smaller or consistent with one in the range of parameters studied, 
with a very weak tendency to increase with increasing lattice size. Due to the ambiguities of 
the estimation of r in the situation of almost complete decorrelation, we do not want to give 
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an estimate for z here. We only state that the results indicate that CSD is almost completely 
eliminated by the time slice blocking algorithm. 

Table 6: Static observables in the two dimensional SU{2) lattice gauge theory on LxL lattices, 

L/e ^ 10. 



L 


statistics 


discarded 


^'^ V16' 16'' 


w^(it) 


P 


p2 


16 


100 000 


10 000 


0.65816(10) 


0.18751(14) 


0.0002(4) 


0.01564(7) 


32 


100 000 


10 000 


0.67917(4) 


0.21287(12) 


0.0002(3) 


0.00855(4) 


64 


50 000 


10 000 


0.68525(6) 


0.2204(2) 


-0.0003(5) 


0.00610(5) 


128 


40 000 


5 000 


0.68688(6) 


0.2222(3) 


0.0008(9) 


0.00545(8) 


256 


40 000 


5 000 


0.68720(7) 


0.2232(2) 


0.0001(6) 


0.00519(4) 



Table 7: Integrated autocorrelation times Tint for the two dimensional SU{2) lattice gauge 
theory on LxL lattices, L/^ ^ 10. If no value is given, we found almost complete decorrelation. 



L 


statistics 


discarded 


^*nt,PF(iA) 


^mt,PF(-|,f) 


Tint,P 


Tint,P^ 


16 


100 000 


10 000 


0.54(1) 


- 


- 


- 


32 


100 000 


10 000 


- 


0.60(1) 


- 


- 


64 


50 000 


10 000 


0.67(1) 


0.70(1) 


0.71(1) 


0.59(1) 


128 


40 000 


5 000 


0.76(2) 


0.74(2) 


0.92(2) 


- 


256 


40 000 


5 000 


0.88(3) 


0.83(2) 


1.01(3) 


- 



9.3 Suminary 

In this section we reported on a multigrid Monte Carlo simulation of SU{2) lattice gauge theory 
in two dimensions. The time slice blocking algorithm with piecewise linear interpolation was 
implemented in the unigrid style with a V-cycle. We used a heat bath version in randomly 
chosen f/(l) subgroups. The simulations were performed on physically large lattices of the size 

m X loe 

All observed integrated autocorrelation times are found to be smaller or consistent with one 
on lattice sizes up to 256^, with a very weak tendency to increase with increasing lattice size. 
Therefore we conclude that the time slice blocking algorithm eliminates CSD almost completely. 

It is fair to say that we used a special feature of two dimensional lattice gauge theory: In 
the infinite volume limit or for open boundary conditions one can decouple two dimensional 
lattice gauge theory to a set of independent one dimensional spin models by choosing the axial 
gauge. Although we use periodic boundary conditions here, our method is very much in the 
spirit of updating on independent time slices. 

It will be discussed in the next section whether this concept is still successful if we generalize 
it to four dimensions. 
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Figure 15: Autocorrelation functions p{t) in the two dimensional SU{2) pure 
lattice gauge theory, j3 = 1024 on a 256^ lattice. Top: p(t) for P, bottom: p(t) 
for P2. 
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10 Multigrid methods for lattice gauge fields in four di- 
mensions 

In this section we generalize the multigrid procedure for pure lattice gauge theory from two 
dimensions to four dimensions. We start with a discussion of the abelian case. 

For nonabelian gauge fields in four dimensions additional difficulties arise because of the 
more complicated geometrical structure of nonlocal updates in four dimensions compared to 
two dimensions. We study the behavior of acceptance rates with increasing block size Lb and 
analyze the special features of the algorithm in the four dimensional case in detail. 

10.1 The abelian case 

The generalization of the multigrid algorithm for U{1) gauge fields of Laursen, Smit and Vink 
(as discussed in section P) from two dimensions to higher dimensions is straightforward: 

Nonlocal updates are formulated as follows: One chooses a hypercubic block x'^ of size L^ 
and a direction r with 1 <t < d. All the link variables U^^t attached to sites x inside the block 
x'g are proposed to be changed simultaneously: 

[4,T -* exp(zs?/'a;)t/x,T , (10.1) 

or in terms of the link angles 

e^,r -^ 0x,r + S^Ij^ . (10.2) 

The kernel ip obeys the normalization condition ( p.8|) . The canonical dimension of the kernel 
ip and of the angle 6 is that of a vector potential, i.e. (2 — d)/2 in d dimensions. The simplest 
choice for ip is again a piecewise constant kernel, 

^const = \ ^B lOr X fc X .^Q g. 

^ [0 for X ^ x'o . 

Let us now study the acceptance rates for these update proposals with the help of formula 
(|3.8|) . We consider general kernels if). For hi = (A7i) we find (cf. appendix^ 

hi = PPT.T.1^- cos(s(^,+^ - V-.))] , (10.4) 

with P = (Trf/p). 

Because of the statistical independence of adjacent time slices we choose ipx to be constant 
in the r-direction, with the support of ip in r-direction restricted to the block x'^. Then we 
obtain 

hi = PPLb E E [1 - cos(s(^,+^ - ^,))] . (10.5) 

For small s, this can be approximated by 

hi ^ \s^(3PLb E E(^-+A - ^x)' = \s^(^Pad-i , (10.6) 



with Ud-i = {tp',—Ailj'). As in two dimensions ip' denotes the kernel ip, restricted to the 

i.[. We absorbed a factor L^^ 



1 /9 

d — 1 dimensional time slice AT. We absorbed a factor L^ in ip' (now normalized as a c? — 1 



dimensional kernel). 
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L, 



Figure 16: Illustration of the geometry of time slice blocking in three dimensional lattice gauge 
theory. 



In the special case of piecewise constant kernels we find for hi 

h = 2{d - 1)/3PL^-^[1 - cos(sLg"'')/')] . 



(10.7) 



From the kinematical point of view, the behavior of acceptance rates in the U{1) lattice gauge 
theory in d dimensions is the same as in massless free field theory. 

This multigrid algorithm with constant interpolation and a W-cycle was studied for the 
simulation of four dimensional U{1) theory by Laursen and Vink |2^. In this model there is 
a deconfinement phase transition between a confinement phase and a Coulomb phase. The 
transition is believed to be of weakly first order |B3[ . The authors report on an acceleration of 



the Monte Carlo dynamics in both phases close to the transition point. This can be understood 
by our kinematical analysis. However the very low frequence of changes from one phase to the 
other on larger lattices in the simulation with a local algorithm could not be accelerated by the 
nonlocal updating. Like the topological problems in the Xy-model and in the two dimensional 
[/(I) model, the problem of metastability can not be resolved from a purely kinematical point 
of view. 



10.2 The nonabelian case: gauge group SU{2) 

10.2.1 Covariant time slice blocking for SU{2) in four dimensions 

We now discuss a generalization of the covariant time slice blocking algorithm of section ^ from 
two dimensions to four dimensions. Many steps of the two dimensional method can be translated 
directly to the four dimensional case. Other features such as the nontrivial background field in 
the bottom of higher dimensional blocks will require a refined treatment. 

Nonlocal updates can be defined as shown in figure ^ (for simplicity illustrated in three 
dimensions). Let us consider a fixed time direction r with 1 < r < 4 and a three dimensional 



time slice AT 



{x G Aq I XT- 
ink vari 
changed simultaneously: 



in AJ". All the link variables U^^t 



t}. One chooses a cubic block x'^ of size L^ that is contained 
attached to sites x inside the block x'^ are proposed to be 



U^,r^Kr = M9)U:. 



(10.^ 
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where Rx{g) = Q^RxQx and g^ G SU{2). The rotation matrices R^ G SU{2) are again 
parametrized as 

Rx{n,s) = cos(s'?/'a;/2) + 2sin(s-?/'2:/2) n-a, (10.9) 

where n denotes a three-dimensional real unit vector, and the ctj are Pauli matrices. ip will 
have support on the three dimensional block x'^. 

Up to now the (7-matrices are arbitrary. In the two dimensional case we have chosen them 
according to the block axial gauge (|8.18|) : Choose the gf- matrices in the one dimensional bottom 
of the block x'^ such that 

Ul^, = gxUx,p,gl+^ = 1 for all (x, x + p,) ex'^ . (10.10) 

The gauge transformed link variables f/J' in the bottom of the block were equal to unity and 
for the case of piecewise constant interpolation the plaquettes in the interior of the block stayed 
unchanged (recall the discussion following figure |Tip. Therefore only the two plaquettes at the 
boundary of a block contributed to the energy change of a piecewise constant update. The 
energy cost was proportional to the surface of the block, not proportional to the volume of the 
block. 

By choosing the block axial gauge in two dimensions we used the fact that we always can 
gauge link variables along a one dimensional open line to unity. Any nontrivial content of the 
gauge field in the one dimensional bottom of the block could be shifted outside of the block by 
applying the block axial gauge. 

However in more than two dimensions, the bottom of a block will contain closed loops (see 



figure |l^). Since a parallel transporter along a closed loop is gauge invariant, we can not get rid 
of the nontrivial curvature that is contained in the loop by any gauge transformation. Therefore 
for nontrivial gauge fields we can not find a gauge transformation g such that U^ = 1 for all link 
variables in the bottom of the block x'^. This means that for nontrivial gauge fields all timelike 
plaquettes that share a link with the bottom of the block will contribute to the average energy 
change of the update. 

The intuitive Leitmotiv for the development of new multigrid methods from section ^ was: 
A piecewise constant update of a nonlocal domain should have energy costs proportional to the 
surface of the domain, but not energy costs proportional to the volume of the domain. 

Unfortunately, according to the discussion above, the nontrivial background field in the 
bottom of the block will lead to energy costs proportional to the volume of the block. Let us 
nevertheless attempt to have as little energy costs as possible: 

Consider the extreme case oi P ^ 00. Then the allowed configurations are pure gauges, i.e. 
configurations that are gauge equivalent to U^^^ = 1 for all x, /i. If we choose g as the trans- 
formation that brings all links to unity, it is obvious that the plaquettes in the interior of the 
block will not be changed by a piecewise constant update. In particular, to have this property, 
it is sufficient to gauge all links inside the bottom of the block to unity. This consideration 
leads to the following gauge condition: Choose g as the gauge transformation that maximizes 
the functional 

Gc,x'^iU,g)= Y. TrigxUx,,g:^^) . (10.11) 

(x,x+fi)Gx'^ 

We call this gauge "block Coulomb gauge" .Q For finite f3 this gauge will not bring all the links 
in the bottom of the block to unity, but still as close to unity as possible. Therefore the gauge 

^ Gauge conditions on the lattice and their relation to transversality conditions for the gauge potential A 
are discussed in appendix O 
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field in the bottom of tlie block is as smooth as possible. This leads to a kind of minimization of 
the energy costs from the interior of the block. Note that in two dimensions the block Coulomb 
gauge condition reduces to the block axial gauge ( |8.18| ). 

The main difference compared to the two dimensional case is that we have to expect that 
the energy change of the update will be proportional to the volume of the block x'^. This 
property is caused by the fact that in four dimensions the gauge field in the bottom of the 
block is smooth but nonzero. This will lead to an algorithmic mass term that grows quadratic 
with the block dimension L^- We are going to investigate the behavior of this algorithmic mass 
in detail below. 

Let us summarize the steps of the nonlocal updating scheme for SU(2) in four dimensions: 

1. Choose a block x'^ of size L^ that is contained in the time slice A[. All link variables Ux,t 
pointing from sites x inside the block in r-direction will be moved simultaneously. 

2. Find the (yf-matrices defined by the block Coulomb gauge condition 

Gc,x'oiU,g)= J2 Tr(5f^[/^,^5f*.^^) = maximal. (10.12) 

(x,x+fj,)£x'g 

3. Propose new link variables U'^^^ by 

U,,r^U',^, = R,{g)U,,r, (10.13) 

with R^{g) = g*R^g^ and 

Rx{n,s) = cos(s?/'a;/2) + isin(s^2:/2) ri-cT. (10.14) 

s is a uniformly distributed random number from the interval [— £:,e], and n is a vector 
selected randomly from the three dimensional unit sphere. 

4. Calculate the associated change of the Hamiltonian A7i and accept the proposed link 
variables with probability min[l,exp(— A7i)]. 

The argument that the detailed balance condition is fulfilled by this updating scheme is 
analogous to two dimensions (cf. section §). The only additional point is that although we only 
have an iterative gauge fixing algorithm in four dimensions, we do not have to fix the gauge 
perfectly. If we always use the same procedure in finding g (e.g. a given number of relaxation 
sweeps starting from (7 = 1), we will always get the same g and the nonlocal update is reversible. 

The detailed implementation and simulation of nonlocal block updates will be described in 



section 11. 



10.2.2 Acceptance analysis for nonlocal S'f/(2)-updates 

First numerical studies revealed that there is no substantial difference in the acceptance rates 
when instead of using the block Coulomb gauge condition one uses the Coulomb gauge condition 
for the whole slice Aj": 

GciU,g)= Yl Tr(5f^f/^,^5f*^^) = maximal. (10.15) 

{x^x+fj^eh-l 
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For a detailed discussion of Landau and Coulomb gauges on the lattice see appendix 0. From a 
practical point of view the Coulomb gauge condition is very convenient: The ^f-matrices can be 
calculated once and then be used for all block sizes Lb- In the block Coulomb gauge they would 
have to be recalculated for every individual block lattice. In addition, the relaxation algorithm 
to determine the ^f-matrices according to the Coulomb gauge condition can be vectorized in a 
straightforward way. 

The energy change associated with the update proposal (|1U.8D is 



A^ = -f ETr(f/^-f/7^) = -f E ETr{(i?:f/f,,i?.+A-f^f,.)^f,;}, (10.16) 

with H*^^ = Ux+fi^TU*_^_^^^U*^^ and U^^^ = g^Ux^^glj^jx- H^ is defined analogously. The relevant 
quantity for the acceptance rates is hi = {AH). For piecewise constant kernels and the gauge 
condition ( |10.15| ) we get (cf. appendix ^ 

hi = 3(3 A {Lb - 1)L| sm\sLB^^y2) + 6^PL|[1 - cos{sLb^^'^ /2)] , (10.17) 



with 



A = {\TT{{Ul^-n.aUl^n.a)Hi:^)) , (10.18) 

P = {^TtUv) . (10.19) 



To the first term in eq. ( |10.17|) all links contribute that are entirely inside the block, whereas 



the second term contains the contributions of all links that have one site in common with the 
surface of the block. For small s, the first term behaves like s^L^. This is exactly the behavior 
of a mass term that, as we have learned in the previous sections, can be toxic for the multigrid 
algorithm. Note that the Coulomb gauge attempts to minimize this mass term by minimizing 
the quantity A. We identify the square root of [3 A with a "disorder mass" mu. 



mD = \J(3A. (10.20) 

To have a physical interpretation of m^ let us discuss the "disorder scale" Id that is given by 
the inverse of the disorder mass: 

mr) \JpA 

If the block size Lb gets of the order of the disorder scale Id, the mass term in eq. ( |10.17| ) 
becomes of order one, and the amplitudes of the nonlocal moves become suppressed. Now the 
crucial question is: How does the scale Id behave for large (3 in comparison with the physical 
correlation length ^? (^ is given by the inverse glueball mass or the inverse square root of the 
string tension.) If Id scaled with ^ the algorithm could efficiently create fluctuations up to the 
scale of ^, as required by the physics of the model. Everything would be all right if for large ^ 

-J- ^^ const . (10.22) 



^oo 



In 4-dimensional SU{2) lattice gauge theory the correlation length is known to increase expo- 
nentially fast with (3. Therefore, for eq. ( [10.22| ) to hold, we would need that the disorder mass 
ttld decreased exponentially fast with increasing (3. Formulated differently, the crucial question 
is whether or not rriD scales like a physical mass. 
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10.2.3 Monte Carlo study of m^ 

We computed mo by Monte Carlo simulations for several values of (3. To maximize Gc we used 
a standard Gauss-Seidel relaxation algorithm vectorized over a checkerboard structure. The 
relaxation procedure consists in going through the lattice and minimizing the gauge functional 
(|10.15|) locally. For production runs it would be advantageous to use an accelerated gauge fixing 
algorithm such as overrelaxation or multigrid |53, p5[. In the Monte Carlo studies reported in 



this section, we always used 50 GauB-Seidel sweeps to determine g. Note that by this procedure 
Gc is not entirely maximized, especially on very large lattices where the relaxation algorithm 
suffers from CSD. However, for the detailed balance to be fulfilled, we only need that one uses 
always the same number of relaxation sweeps. Several tests revealed that increasing the number 
of relaxation sweeps beyond 50 did not affect the acceptance rates in a substantial way. In our 
implementation, 50 Gaufi-Seidel sweeps over all slices of a given direction r required the same 
CPU time (on a CRAY Y-MP) as four Creutz heat bath SU{2) update sweeps. 



n(s 



1.0 



0.5 



0.0 




0.0 



0.5 



1.0 



Figure 17: Q{s) in 4-dimensional SU{2) lattice gauge theory using piecewise 
constant kernels, (3 = 2.6 on a 20^-lattice. From top to bottom: Lb = 2,4,8, 16. 
Points with error bars: Monte Carlo results, lines: analytical results using m,£, 
and P from Monte Carlo (errors smaller than line width). 



We checked the validity of the acceptance formula (|3.8| ) using Monte Carlo estimates for 
niD and P. Figure |1^ shows results for P = 2.6 on a 20"^ lattice. The results perfectly justify 
the usage of the approximation formula. It is therefore sufficient to study the behavior of the 
quantities mo and P. Our Monte Carlo results are presented in table ^ The last column 
gives the statistics in sweeps (equilibration sweeps are not counted here). We used a mixture 
of four microcanonical overrelaxation sweeps followed by a single Creutz heat bath sweep. 
Measurements (including the determination of g) were performed every 25 sweeps. 
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Table 8: Monte Carlo results for mj:, and P 



lattice size 


P 


mo 


P 


statistics 


8"^ 


2.4 


0.507(2) 


0.6305(3) 


10 000 


12^ 


2.4 


0.4957(4) 


0.6300(2) 


10 000 


16'* 


2.4 


0.4955(2) 


0.62996(5) 


10 000 


8* 


2.6 


0.497(4) 


0.6703(1) 


30 000 


W 


2.6 


0.465(2) 


0.6702(1) 


20 000 


IQ^ 


2.6 


0.4644(3) 


0.67004(5) 


10 000 


20'* 


2.6 


0.4650(2) 


0.67008(5) 


5 000 



Table 9: Comparison of rriD with physical masses 



lattice size 


P 


TUd 


y/K 


mo+ 


mo/y/K 


mD/mo+ 


16* 


2.4 


0.4955(2) 


0.258(2) 


0.94(3) 


1.92 


0.53 


20* 


2.6 


0.4650(2) 


0.125(4) 


0.52(3) 


3.72 


0.89 



In table |^ we display the ratios of the disorder mass mo with two physical masses, the square 
root of the string tension k and the lowest glue ball mass mo+ . The estimates for the physical 
masses are taken from ref. |EB|. The results show that the disorder mass is nearly independent 
of P in the range studied, whereas the physical masses decrease by roughly a factor of two. 
Thus, rriD is not scaling like a physical mass for the couplings studied here. We conclude from 
this that for large blocks the term quadratic in Lb will strongly suppress the acceptance rates 
even when the ratio of correlation length and block size Lb is kept constant. 

From this kinematical analysis it is clear that we can not expect that CSD will be eliminated 
by such nonlocal updates. 

Let us give a plausibility argument for the large f3 behavior of m.^ = y/JJA: From the 
definition ( |10.18| ) of A it is clear that A vanishes for large /? because U^^^ goes to unity in this 
limit. Since A is a quantity that is defined on the scale of the plaquette, it is dominated by 
the local disorder. A has nothing to do with collective excitations of the gauge field that are 
responsible for the formation of glueballs with a mass that decreases exponentially in f3. The 



leading weak coupling behavior of the plaquette is |E7| 



P 



^-ii^^^fc 



Therefore it is natural to expect a weak coupling behavior for A like 

A = 

with a constant c. If this conjecture was true we would get 

m,£) = JpA — > const . 




(10.23) 



(10.24) 



10.25) 
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In four dimensions the range of j3 values and lattice sizes studied here is too small to check the 
large /3 behavior of m/5 in detail. Let us study the analogous situation in two dimensions. 

10.2.4 Digression: Monte Carlo study of mo in two dimensions 

We want to study the kinematical effect of a remaining nontrivial background field U^^^ for 
/i 7^ r in the bottom of a block in two dimensions. This is of course an artificial complication 
of the situation in two dimensions. We always can choose the block coulomb gauge there, as 
we did in the previous sections. Nevertheless we can use the two dimensional case as a simple 
example where the region of large (3 is easily accessible, in contrast to four dimensions. In 
analogy to the four dimensional case we are going to take the (yf-matrices within a time slice 
according to the Coulomb gauge condition 

GciU,g)= Y^ Tr(5f^f/^,^5f*+^) = maximal . (10.26) 

(x,x+/i)GA[ 

In two dimensions this maximization has a simple solution: Take the ^f-matrices such that 

Ul, = 9.U,,,9Uf. = Pl^' for all x G A[ . (10.27) 

Pt is the Polyakov loop in space direction, 

L 

Pt= U U.^f. (10.28) 

on a L X L lattice, and P^ is the closest root to unity. 

By taking the Coulomb gauge condition in two dimensions we now have a smooth but 
nontrivial background field in the bottom of the block. Let us study the effect on the kinematics. 
For piecewise constant interpolation we have (cf. appendix ^ 



1/2 



1/2 



hi = 2(3P [I - cos [sL'jl^")\ + I3A{Lb - 1) sin' (sL^^'j , (10.29) 

with 

A = (iTr((t/f,^ - n-aUl^n-aW:,)) , (10.30) 

where g is taken according to eq. ( |10.26| ). 

We computed A and mo = \/l3A by Monte Carlo simulations in the region of large (3 



close to the continuum. Our Monte Carlo results are presented in table ^ The statistics are 
counted in sweeps (we started from equilibrated configurations). We used a mixture of four 
microcanonical overrelaxation sweeps followed by a single heat bath sweep. Measurements were 
performed every 5 sweeps. 



The large (3 behavior of A is shown in figure 0. It is consistent with A ~ 1/13. In table 
pn] we display the ratios of the disorder mass mo with the physical mass, the square root of 
the string tension k. The values for the string tension are taken from the exact solution in the 
infinite volume limit (appendix P). The results show that the disorder mass in two dimension 
goes to a constant for large /?, whereas the physical mass decreases hke 1/a/5. Therefore, rriD 
is definitely not scaling like a physical mass in the continuum limit in two dimensions. 
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0.01 



0.008 



0.006 



0.004 



0.002 




Figure 18: Dependence of the quantity A on 1/(3 in two dimensional SU{2) 
lattice gauge theory. Diamonds: Monte Carlo results. The line corresponds to 

A ~ 1//3. 
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Table 10: Monte Carlo results for P, A and mn in two dimensions 



lattice size 


P 


P 


A 


mu 


statistics 


82 


1 


0.2400(4) 


0.01357(15) 


0.1165(7) 


20 000 


162 


4 


0.65791(14) 


0.00951(3) 


0.1950(4) 


20 000 


322 


16 


0.90783(2) 


0.003284(6) 


0.2292(2) 


20 000 


642 


64 


0.976655(3) 


0.000887(1) 


0.2383(2) 


20 000 


1282 


256 


0.9941462(3) 


0.000225(2) 


0.2400(1) 


20 000 


Table 11: 


Comparison of mn with physical masses in two dimensions 




lattice size 


P 


e 


TTlD 


^Jk 


■niD/s/K. 






82 


1 


0.8373 


0.1165(7) 


1.1942 


0.0976 






162 


4 


1.5458 


0.1950(4) 


0.6469 


0.301 






322 


16 


3.2155 


0.2292(2) 


0.3110 


0.737 






642 


64 


6.5065 


0.2383(2) 


0.1536 


1.55 






1282 


256 


13.0511 


0.2400(1) 


0.0766 


3.13 





We conclude from this that for large blocks the term quadratic in Lb will strongly suppress 
the acceptance rates also in two dimensions. Therefore the smooth but nontrivial background 
field in the coulomb gauge that is still left in the one dimensional bottom of the time slice will 
cause CSD in this variant of the algorithm. 

In analogy to this result in two dimensions we expect the same behavior of niD in the large 
j3 limit in four dimensions. 

10.2.5 Comparison of volume and surface effects in four dimensions 

By the kinematical analysis in four dimensions we found that CSD has to be expected because 
of a mass m^i generated by the local disorder in the bottom of the blocks. However, one could 
hope that the value of the unwanted mass term is so small that it is only harmful at very large 
values of (3 on huge lattices. Let us examine the effect of this term in more detail. 

Recall that hi is built up from two contributions. The first contribution is that related to 
the gauge field disorder inside the block and is quantitatively represented by the mass mo- 
The second contribution is associated with the block surface. The latter can of course be made 
smaller by using smooth kernels tp instead of the piecewise constant kernels discussed so far. 
However, the disorder term cannot be expected to become smaller for smooth kernels (see 



below). In figure |T9| we plotted separately the two contributions to hi 



hi,A 
hiP 



3(3A {Lb - l)Ll sin^isL's^^'^ /2) 
Q(3P Ll{l - cos{sLb^^^ /2)) , 



;io.3i) 



for j3 = 2.6 and block size Lb = 8 on a 20'^ lattice. The plot shows that already for this block 
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size the disorder contribution is by no means small - it is comparable to the surface effect. It is 
therefore not clear that one could achieve any significant improvement by using smooth kernels. 
To investigate this in more detail, we derive an expression for hi (cf. appendix ^, valid for 
smooth kernels as well: 






[10.32) 



Since I]^x ~ ^b^ ^^ S^^ essentially the same behavior for the disorder contribution as in the 
case of piecewise constant kernels. 



4.0 



hi,A(s), hi,p(s) 



2.0 



0.0 




0.0 



0.2 



0.4 



0.6 



Figure 19: Comparison of disorder and surface effects for four dimensional SU{2) 
lattice gauge theory using piecewise constant kernels on an 8^-block, (3 = 2.6 on 
a 20'^-lattice. Solid line: hi^A^s) (disorder effects), dashed line: /ii,p(s) (surface 
effects). 



We checked this quadratic approximation against the exact expression ( |10.17|) in the case 
of piecewise constant kernels. For Lb = 8 the deviations are already negligible in the range of 



s-values displayed in figure ^ For smooth -0"^ 
two contributions to hi 



kernels we show separately in figure 20 the 






hi,p.A = ^s^{P-A)a3. 



10.33) 



for f3 = 2.6 and block size Lb = 8 on a 20^ lattice. We observe that the surface effects are 
lowered by the smooth kernels, but the disorder contribution is even higher than for piecewise 
constant kernels. 
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J-Ila^^/) J-1i,]p-a\^ 
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Figure 20: Comparison of disorder and surface effects for the four dimensional 
SU{2) lattice gauge theory using smooth ■?/;**"<^ kernels on a 8^-block, (3 = 2.6 on 
a20^-lattice, quadratic approximation used. Solid line: hi^A^s) (disorder effects), 
dashed line: hi^p^A{s) (surface effects). 



Piecewise constant kernels have the practical feature that once the change of the Hamil- 
tonian has been calculated, one can perform multihit Metropolis updating or microcanonical 
overrelaxation. In a special case, even a recursive multigrid implementation with a W-cycle is 
possible (see section |ll] below). For smooth kernels the change in the Hamiltonian would have 
to be calculated again and again. Also the advantages of smooth kernels are not that clear 
on small three or four dimensional blocks. For the concrete simulation we will use piecewise 
constant kernels. 

10.2.6 Maximally abelian gauge 

Our proposal for the choice for g was motivated by the desire to minimize the quantity A. We 
now ask whether there is a better choice for g than the g determined by the Coulomb gauge 
condition. For the sake of simplicity let us take n = (0, 0, 1), i.e. n-a = a^. Then A is given by 

A = (iTr((f/f,^ - as f/f,, as)//^;)) . (10.34) 

The choice of the Coulomb gauge condition aimed at bringing f/|^ as close to unity as possible. 
Alternatively, one might require that f/|^^ should be as close as possible to a S'f/(2)-matrix of 
the form oq + icLsa^. This will also lead to a small A. The corresponding gauge transformation 
g can be found by maximizing the functional 

GAiU,g)= J2 Tr(a3f/,V3^,^;) , (10.35) 
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leading to the inaxiinally abelian gauge |5^, here implemented only on a slice. We computed 



rriD also using the g^s resulting from this gauge condition and compared the results with the 
ones obtained by using the Coulomb gauge condition. We did not find a substantial difference. 
We prefer the Coulomb gauge condition because it does not depend on the direction n and thus 
saves computer time. 

10.3 Summary 

In this section we discussed the generalization of multigrid Monte Carlo methods for the sim- 
ulation of pure gauge fields from two dimensions to four dimensions. 

In the abelian case the generalization is straightforward. We have seen that the kinematical 
behavior of the resulting method is similar to massless free field theory. Therefore a considerable 
acceleration can be expected and was indeed observed in four dimensional U{1) lattice gauge 
theory [^. However, metastability causes additional problems close to the phase transition 
which can not be understood in terms of our method. 

In the second part of this section we have generalized the time slice blocking algorithm that 
lead to successful acceleration of the simulation of SU{2) gauge fields in two dimensions to four 
dimensions. In four dimensions, new difficulties occur that are due to the nontrivial background 
field in the bottom of the blocks. We investigated the scale dependence of acceptance rates 
in detail. Here we found that an algorithmic mass term generated by the disorder in the 
background field suppresses the acceptance rates on large blocks. From our kinematical analysis 
we can not expect that the proposed algorithm will have a chance to reduce the dynamical 
critical exponent below z k.2. However, compared to local Monte Carlo algorithms there could 
still be an acceleration of the dynamics by a constant factor, depending on the details of the 
implementation. This question will be investigated in the following section. 
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11 Multigrid Monte Carlo simulation of SU{2) lattice 
gauge fields in four dimensions 

In this section we study whether our nonlocal Monte Carlo updating can accelerate the Monte 
Carlo dynamics of a local heat bath algorithm by a constant factor. We do not aim at estimating 
the dynamical critical exponent z. 

We begin with a description of the implementation of the time slice blocking scheme for 
SU{2) in four dimensions as a recursive multigrid algorithm with piecewise constant interpola- 
tion and a W-cycle. Then we give a detailed report on a comparison of the time slice blocking 
algorithm with the local heat bath algorithm. 

11.1 Implementation of the time slice blocking 

For a concrete implementation of the multigrid Monte Carlo algorithm in four dimensions we 
choose the time slice blocking method: In a basic time slice blocking step, all links pointing 
from a given three dimensional time slice AJ" in the r-direction are going to be updated. First 
we are going to describe the implementation of such a basic step for a given time slice A[. In 
particular we explain a recursive multigrid procedure that can be constructed by updating in 
a global f/(l) subgroup of SU{2). How time slice blocking steps are applied to different time 
slices A[ and to different orientations r in sequence will be explained in a second paragraph. 

11.1.1 Organization of a basic time slice blocking step 

We now describe the time slice blocking method for a fixed three dimensional time slice AJ" = 

{a; G Ao I Xt- = t}. 

The basic updating step is based on a local, site dependent "rotation" of all link variables 
pointing from sites s G A[ in the r-direction by 

U,,r^U'^^^ = R,{g)U,,r, (11.1) 

with R^{g) = glRxQx and 

Rx{nx,Ox) = cos{6x) +ism{9x)nx-a. (11.2) 

The gauge transformation g is obtained by imposing the Coulomb gauge condition ( |10.15|) on 



slices A[ as defined in section IT^ 



We are interested in a recursive multigrid implementation. Therefore we use piecewise 
constant kernels and update a fixed U{1) subgroup of SU{2) globally. In this special case an 
explicit multigrid scheme can be constructed. This feature was explained in ref. [^ in a related 
context. 

The restriction to a global U{1) subgroup is defined by taking the directions fix of the em- 
bedded f/(l)-subgroups to be independent of the site x within the time slice, i.e. fix = n for 
all X G A[. The global n direction is taken as a random vector from the three dimensional 
unit sphere. Then a conditional Hamiltonian on the finest lattice can be calculated by substi- 
tuting the "rotated" gauge field U' in the fundamental Hamiltonian. The resulting conditional 
Hamiltonian 7io{S) is of the type of a generalized three dimensional XF-model: 

- HoiS) = E E S^Bx,^Sx+f. + const , (11.3) 
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where S^ denotes the Xy-spin variable 



«- 3: ■ (^i'^) 



and S]^ its transpose. B^^fj, is a hnk dependent real 2x2 matrix 



x,ii X,fl 



Bx,fj. — I j^21 j^22 ) (11-5) 



with the matrix elements 



BZ = fTr(f/f,^^n-afff;) , 

Bll = -|Tr(m-af/f,^if,^;) , 

BZ = -^TT{zn.aUl^in-aHi:;) . (11.6) 

Now we set up a recursive multigrid procedure to update the three dimensional generalized 
XF-model on a time slice A[. The first blocking step is performed by dividing the fine lattice A[ 
in cubic blocks x' of size 2^. This defines a three dimensional block lattice Ai. Then a nonlocal 
block update consists of a constant rotation of all XF-spins S^ within a three dimensional 
block x' G Ai by the block angle 6x''- 

An equivalent parametrization of this update is 

Sx ^' S^ = Aj.Sx' , (11-8) 

where the 5'0(2)-matrix A^ and the block-XF-spin S^' are given by 

_ / cos(^^) -sin(^^) \ _ / cos(^^/) \ , . 

^^ - 1^ sin(^J cos(^J j ' ""^'-y sin(^.0 ) ' ^'^'^^ 

By iterating this procedure one gets conditional Hamiltonians 7ik{S) on three dimensional 
coarser layers A^. An important point is that in the special case considered here 7ik{S) always 
stays of the form 

- HkiS) = E E S^B^^^S^^^ + E S^M^S, + const (11.10) 

with 2x2 coupling matrices B^^^ and M^ that depend on the links {x, ^) and on the sites x in 
Afc respectively. The coupling matrices B^'^n and M^' that depend on the block links {x',fi) or 
blocks x' on the next coarser lattice A^+i can be calculated recursively: 

Bx',fj, = 2^ A^B^^^Ax-^-fi 

X G x' 
x+fiGx'+(i 

M^, = E AlBx^^Ax^f,+ Y.AlMxA,. (11.11) 

xex' x^^' 
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Here x' + fi denotes the next neighbor block of the block x' in the coarse layer Afc+i. Note 
that by the recursion the conditional Hamiltonian 7ik{S) on a coarser layer A^ contains not 
only next neighbor couplings of the XF- variables (as the conditional Hamiltonian 7io{S) in 
eq. ( |11.3|) on the finest lattice) but also site-site couplings of the type of a mass term with site 
dependent mass matrix M^. The algorithmic mass term discussed in section ^ that leads to 
a decrease of acceptance rates on large blocks originates exactly from this site-site term in the 
conditional Hamiltonian. 

This recursive blocking procedure is repeated until we end up with the coarsest layer that 
consists of a single point. (I. e. in our simulation on a 8^-lattice we are going to update on 
block lattices of the size 4^, 2^ and 1^.) 

After a blocking step to a coarser layer Ak+i the initial configuration on this coarse layer is 



e^, = ^ S^> = { ^ ] for all x' E A^+i . (11.12) 




On each layer we perform one sweep of 10 Hit Metropolis updates in terms of the 6^ variables: 

e,-^e'^ = e, + Ae , (11.13) 

where A9 is taken to be equally distributed from the interval [—£,£]. We tuned the Metropolis 
step size e in order to obtain acceptance rates of ~ 50% on each layer. If we return from a 
coarse layer A^+i to the next finer layer A^, the configuration {S} on the finer layer is updated 
according to eq. ( |11.7[ ). 

In this implementation, the sequence of visiting the different three dimensional layers of a 
time slice is organized according to a W-cycle. From the point of view of numerical work this is 
possible because of the recursive definition of the algorithm. In principle the described imple- 
mentation is not only feasible in the time slice blocking formulation but also with hypercubic 
four dimensional blocks. 

11.1.2 Sequence of basic time slice blocking steps 

The sequence of the updating on different time slices is as follows: On an L^-lattice we visit 
the L time slices Aj, t = 1, . . . L. In this way, all link variables Ux,i on the lattice that point 
in the 1-direction are updated in a nonlocal way on block lattices with Lb = 2,4, . . .. Then 
we perform a local Creutz heat bath sweep through all links on the lattice. Now we change 
the time direction r from r = 1 to r = 2: We visit the L time slices A^, t = 1,. . .L, again 
followed by a sweep of local heat bath updates. The same scheme (a visit of all time slices and 
a local heat bath sweep) is repeated also for the r = 3 and r = 4 direction, such that all the 
link variables Ux,t have been updated in a nonlocal manner for all r = 1, . . .4. This sequence 
is repeated periodically. 

Observables are measured after each local heat bath sweep. In addition, we perform a 
random translation of the lattice after each local heat bath sweep in order to avoid effects from 
fixed block boundaries [^ . 



To reduce the computer time needed for gauge fixing we use a slight modification: If all 
the link variables Ux,fj, with fi ^ t lying in the bottom of the time slice A[ have been gauged 
according to the Coulomb gauge condition, we update not only the link variables Ux,t pointing 
from this time slice in the positive r-direction, but we update also the links pointing from this 
time slice to the negative r-direction. Thus, only every second time slice A[ has to be gauged. 
This leads to a decrease of the numerical work required for gauge fixing by a factor of two. 
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11.2 Simulation and Results 

11.2.1 Observables 

The observables measured are square Wilson loops 

W{I,I) = {^,Tt{U{Cjj)), (11.14) 

where U{Cij) is the parallel transporter around a rectangular Wilson loop Cjj of size J x /. On 
the SMattice we measure W{1, 1), W{2, 2), and Vr(4, 4). Another important class of quantities 
is built up from timelike Polyakov loops. A Polyakov loop at the three dimensional spatial 
point X is defined by 

L 

P^=iTrn%,t),4. (11.15) 

t=i 

We measure the lattice averaged Polyakov loop 

p=[^j:p^). (11-16) 

the lattice averaged Polyakov loop squared 

and the sign of the lattice averaged Polyakov loop 

sign(P) = /sign l^-l ^ P^^l \ . (11.18) 

11.2.2 Run parameters and static results 

In order to study the acceleration by the multigrid algorithm, we compare the autocorrelations 
of the nonlocal algorithm with a standard local Creutz heat bath algorithm. Precise measure- 
ments of autocorrelation times r require high statistics simulations with run lengths of at least 
lOOOr — 10 OOOr. For reasons of computer time we decided to simulate on relatively small lat- 
tices. The algorithms are compared on a 8^-lattice at /3 = 2.2, 2.4 and 2.6. We started the local 
heat bath runs from ordered configurations (all link variables set equal to unity) and discarded 
a suitable number of iterations for equilibration. For the multigrid simulations we used warm 
starts from already equilibrated configurations. Measurements were taken after each local heat 
bath sweep. 

In order to reduce the computer time needed for gauge fixing, we use 10 sweeps of overrelax- 
ation [Q with overrelaxation parameter uj = 1.7. The degree of maximization of the Coulomb 
gauge functional Gc ( |10.15| ) was found to be about the same as if we used 50 GauB-Seidel 
relaxation sweeps. 

With these run parameters the computer time needed by our implementation on a CRAY 
Y-MP for one measurement on the 8^-lattice by the multigrid procedure is about a factor of 
2.8 larger than the time needed for one measurement by the Creutz heat bath algorithm. This 
factor could still be lowered by using a multigrid method for gauge fixing |]55[] . 



The static results of our runs are given in table |T2|. The values for the heat bath and for 
the multigrid algorithm agree within errors. Where a comparison is possible, our results for 
the Wilson loop observables are consistent with existing data in the literature [^ . 
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Table 12: Comparison of numerical results for various static observables of heat bath (HB) 
and multigrid (MG) simulations on the 8"^ lattice. 



/5 


2.2 


2.4 


2.6 


algorithm 


HB 


MG 


HB 


MG 


HB 


MG 


statistics 


100 000 


50 000 


100 000 


100 000 


100 000 


100 000 


discarded 


10 000 


equi. 


10 000 


equi. 


10 000 


equi. 


WilA) 


0.56922(5) 


0.56917(9) 


0.63024(3) 


0.63026(3) 


0.67029(2) 


0.67029(2) 


W{2,2) 


0.13628(7) 


0.13625(10) 


0.22319(11) 


0.22330(13) 


0.28875(7) 


0.28877(5) 


iy(4,4) 


0.00127(3) 


0.00132(4) 


0.01371(7) 


0.01379(6) 


0.03762(8) 


0.03770(6) 


P 


0.00006(13) 


0.00015(18) 


-0.0020(16) 


-0.00006(13) 


0.008(8) 


0.002(6) 


p2 


0.000584(3) 


0.000580(4) 


0.00211(5) 


0.00202(5) 


0.0092(2) 


0.0096(2) 


sign(P) 


0.0000(2) 


0.0002(3) 


-0.004(3) 


-0.000(2) 


0.013(13) 


0.003(11) 



11.2.3 Autocorrelation times 



Table 13: Comparison of the exponential autocorrelation times of heat bath (HB) and multigrid 
(MG) simulations for different observables on the 8^ lattice. 



/3 


2.2 


2.4 


2.6 


algorithm 


HB 


MG 


HB 


MG 


HB 


MG 


statistics 


100 000 


50 000 


100 000 


100 000 


100 000 


100 000 


discarded 


10 000 


equi. 


10 000 


equi. 


10 000 


equi. 


Tiy(i,i) 


6.9(2.9) 


10.3(3.5) 


5.7(1.5) 


13(9) 


1.8(6) 


1.9(1.0) 


'TW{2,2) 


6.9(1.7) 


8.0(1.8) 


26(6) 


35(7) 


4.0(1.2) 


2.9(6) 


Tiy(4,4) 


^ 1.3 


^ 1.3 


22(3) 


26(3) 


10.4(1.8) 


10.3(2.4) 


-rp 


3.3(5) 


3.8(6) 


93(13) 


67(4) 


279(45) 


274(33) 


Tp2 


1.0(3) 


^ 1.3 


35(4) 


32(6) 


48(8) 


46(6) 


Tsign(P) 


5.3(1.2) 


5.2(1.3) 


92(13) 


68(6) 


275(41) 


277(39) 



The results for the autocorrelation times are given in table |T^ For the observables with 
the longest autocorrelation times P and sign(P), we give a comparison of the autocorrelation 
functions of the multigrid algorithm with the heat bath algorithm for /3 = 2.2, 2.4 and 2.6 in 



figures pTl pg and p3 



Note that the autocorrelation functions p{t) do not in general show a pure exponential 
decay but exhibit a crossover from a fast mode to a slow mode that eventually governs the 
asymptotic decay for large t. Therefore the measurement of the integrated autocorrelation 
time Tint with a self consistent truncation window of 4rj„f might not capture the asymptotic 
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decay of p(t) correctly. In the present case we decided to extract a rough estimate of the 
exponential autocorrelation time Tg^p and the corresponding errors by the following procedure: 
We plotted the p(t) with error bars on a logarithmic scale and decided at what value of t the 
asymptotic exponential behavior started. Then we drew the highest and lowest straight lines 
that were compatible with the data in the asymptotic regime. The errors are given such that 
the highest and the lowest value of r^^p lie at the ends of the interval result ± error. 

All our results for the exponential autocorrelation times Tg^.^ of the multigrid and the heat 
bath algorithm are compatible within errors except from r^xp for P and sign(P) at j3 = 2.4, 
see also figure ^ In this case the superficial gain of the multigrid algorithm is a factor of 
1.4. However the computational overhead of the multigrid method compared to the heat bath 
algorithm (a factor of 2.8 in our implementation) was not yet taken into account. So there is 
no net gain in computer time also in the case of /? = 2.4. 

11.3 Summary 

In this section we studied an implementation of the multigrid Monte Carlo algorithm in four 
dimensions. We used the time slice blocking method and piecewise constant interpolation with 
a W-cycle in a recursive multigrid version by updating in f/(l) subgroups of SU{2). Simulations 
with gauge couplings (3 = 2.2, 2.4 and 2.6 were performed on an 8^-lattice. 

Apart from a modest acceleration of Polyakov loop observables at /3 = 2.4 by a factor of 1.4, 
no improvement was found compared to a local heat bath algorithm. Since the nonlocal update 
procedure has a computational overhead of a factor of 2.8 on the 8^-lattice (on a CRAY Y-MP), 
there is no net gain but even a loss in CPU-time. This factor depends however on the details 
of the implementation. 

Since from the theoretical analysis (section |10|) one can not expect a lower dynamical critical 
exponent z, there is no hope that on larger lattices the method will perform better. 

Possible improvements of our nonlocal updating scheme were investigated by Gutbrod |]BD[ . 



Starting from the Coulomb gauge, he uses an additional smoothing by taking a quadratic ap- 
proximation to the action and updating in terms of approximate eigenfunctions. The resulting 
nonlocal updates are performed as (approximately) microcanonical overrelaxation steps. Asym- 
metric lattices of size L^T with T y> L and anisotropic couplings (e.g. (3 = 4.5 in the time 
direction and /5 = 3.0 in the space directions) are used. The results indicate a superficial gain 
of a factor of 1.5 — 3 compared to a local overrelaxation algorithm. However if the (implemen- 
tation dependent) computational overhead of the nonlocal method is taken into account, the 
net gain is marginal. 
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Figure 21: Comparison of the autocorrelation functions p{t) for the heat bath 
(HB) and multigrid (MG) algorithm for SU(2) on the 8*-lattice at {3 = 2.2. Top: 
Polyakov loop P, bottom: average sign of Polyakov loop sign{P). 
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Figure 22: Comparison of the autocorrelation functions p{t) for the heat bath 
(HB) and multigrid (MG) algorithm for SU(2) on the 8*-lattice at {3 = 2.4. Top: 
Polyakov loop P, bottom: average sign of Polyakov loop sign{P). 
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Figure 23: Comparison of the autocorrelation functions p{t) for the heat bath 
(HB) and multigrid (MG) algorithm for SU(2) on the 8*-lattice at {3 = 2.6. Top: 
Polyakov-loop P, bottom: average sign of Polyakov loop sign{P). 
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12 Summary and Outlook 

We presented a simple yet accurate formula that expresses acceptance rates for nonlocal update 
algorithms in terms of the average energy change hi = (A7Y) of the proposed move. The 
quantity hi depends only on one local observable (or two in four dimensional SU{2) lattice 
gauge theory) that is simple to compute, e.g. by Monte Carlo simulations on a small lattice. 

With the help of our formula we found that all investigated models fall into two kinematical 
classes: For Sine Gordon, 0^ theory and SU{2) lattice gauge theory in four dimensions, s had 
to be rescaled like 1/Lb for piecewise constant and for smooth kernels in order to maintain Lb- 
independent acceptance rates. For massless free field theory, the XY model, the 0{N) nonlinear 
cr-model, the CP^~^ model, U{1) lattice gauge theory and SU{2) lattice gauge theory in two 
dimensions one can achieve s ~ const by choosing smooth kernels. 

The kinematical behavior of multigrid algorithms in interacting models can be related to 
the dynamical critical behavior. By comparing the behavior of the acceptance rates in in- 
teracting models close to a critical point with free field theory, where CSD is known to be 
eliminated by a multigrid algorithm, we found the following rule: For an interacting model, 
sufficiently high acceptance rates for a complete elimination of CSD can only be expected if 
hi = (7i(0 + sip) — 'H{(f))) contains no algorithmic mass term ~ s^J^xi^x- To our knowledge, 
all numerical experiments with multigrid Monte Carlo algorithms that have been performed so 
far are consistent with this rule. With the help of this criterion it is possible to decide whether 
a given statistical model is a natural candidate for multigrid Monte Carlo or not. 

Heuristically, the meaning of the rule can be stated differently: A piecewise constant update 
of a nonlocal domain should only have energy costs proportional to the surface of the domain, 
not energy costs proportional to the volume of the domain. 

To check the prediction that too small acceptance rates on large blocks cause CSD, we 
performed a multigrid Monte Carlo simulation of the Sine-Gordon model in two dimensions. 
We found z = 1.9(1) with piecewise constant interpolation and a W-cycle. It was studied 
whether one can compensate for too small update steps on large blocks by accumulating many 
of these steps randomly. To check this possibility, we simulated the model with a higher cycle 
with 7 = 4. We found z = 1.9(1) also for this variant of the algorithm. 

A multigrid algorithm for nonabelian gauge fields, the time slice blocking, was developed. 
It is based on the statistical decoupling of adjacent time slices as long as only link variables in 
the time direction are updated. The kinematical analysis predicted a strong reduction of CSD 
in a simulation of SU{2) in two dimensions. Almost complete elimination of CSD was indeed 
observed in numerical experiments on lattice sizes up to 256^. 

The multigrid Monte Carlo algorithm was generalized from two dimensions to four dimen- 
sions. In four dimensions, a nontrivial background field in the bottom of the blocks caused 
difficulties. We found that this background field generates an algorithmic mass term that 
suppresses the acceptance rates on large blocks. 

We argued that the algorithm will only work well if the algorithmic disorder mass m^) scales 
like a physical mass, i.e. if rriB) decreases exponentially fast with increasing /5. However, there 
is numerical and theoretical evidence that rriB ~ const for large f3. Therefore the expectation 
is that the proposed algorithm will not have a chance to reduce CSD. 

Whether the nonlocal algorithm is at least able to accelerate the local heat bath algorithm by 
a constant factor was investigated by simulations of SU{2) lattice gauge theory on an 8^-lattice. 
Superficially, a modest acceleration for Polyakov loop observables was found aX (3 = 2.4, however 
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when the additional computational work is taken into account, there is no net gain. From the 
theoretical analysis it can not expected that on larger lattices the method will perform better. 
In summary, it is fair to say that for nonabelian pure gauge theory in four dimensions the local 
overrelaxation algorithm is superior to nonlocal updating schemes up to now. 

The kinematical mechanism that leads to a failure of multigrid algorithms is well described 
by our analysis. As an example, the analysis of the time slice blocking algorithm for nonabelian 
gauge fields in four dimensions can help to understand the additional difficulties one has to face 
in four dimensions compared to two dimensions. We hope that a better understanding can lead 
to improved multigrid algorithms that can overcome kinematical obstructions stemming from 
an algorithmic mass term. 

There are two possible strategies for improved algorithms: 

First one could consider the stated rule as a (non rigorous) no-go theorem. Then one can 
try to avoid kinematical problems beforehand by developing algorithms that do not have energy 
costs proportional to the volume of a block. A natural way is to take a global symmetry of the 
model under consideration and to apply the corresponding symmetry transformation locally on 
blocks. 

From this point of view it is clear that updates that are well adapted to the 0^ model (with 
a discrete Z2 symmetry of the Hamiltonian) are not continuous, but discrete. In fact, the 
very efficient cluster algorithm that has been developed for this model is based on the discrete 
Z2 symmetry |61|. Another example where discrete updates appear to be more natural than 
continuous updates is the Sine Gordon model (with a discrete TL symmetry of the Hamiltonian). 
Generalizations of the cluster algorithm for the discrete Gaussian model |^ (that has the 
same symmetry as the Sine Gordon model) seem to be efficient also in this case [^ . The four 
dimensional nonabelian SU{N) lattice gauge theory has both a global discrete Z^r symmetry 
and a local SU{N) gauge symmetry (which includes a global SU{N) symmetry as a subgroup). 
Which symmetry could be used for a fast nonlocal updating algorithm (if there exists any) is 
an open question. 

A second, more ambitious way to deal with kinematical obstructions is not to avoid them 
but to attempt to overcome them directly. Recently, very efficient methods for the acceleration 
of tunneling processes in systems close to a first order phase transition have been developed ||6^ 



There, the computational effort could be reduced from an exponential to a polynomial growth 
with the volume. 

Since the mass term problem in models with a discrete symmetry can also be interpreted 
as slow tunneling between different minima of a conditional coarse grid Hamiltonian, possible 
combinations of these new developments with multigrid methods might be promising also for 
reducing CSD. However, since CSD is "only" a polynomial problem, the requirement that the 
computational work should be proportional to the volume (and not to a power of the volume 
as in systems close to a first order phase transition) is quite restrictive. 
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A Details of the acceptance analysis for gauge theories 

In this appendix we discuss the details of the calculation of hi = ( A7Y) for gauge theories. 

A.l Abelian gauge fields 

The Wilson action for abelian gauge fields in d dimensions can be written in link angles 6^^^ 

no) = /? E E [1 - cos(V(^))] , (A.l) 

with the plaquette angle 

All the link angles O^^r that point from sites x inside the hypercubic block x'^ in the time 
direction are proposed to be changed simultaneously: 

Ox,r -^ O^^r + Sipx ■ (A. 3) 

The associated change of the hamiltonian is 

AH = -/? E E {cos (V(^) + sii^x+fi -ipx))- cos (V(a;))} 

= ~^ E E {cos(6'^^(a;)) [cos(s(^^.+^ - tp,^)) - 1] - sin(6'^^(a;)) sin(s(^^+/i - ^^.))} • (A.4) 

Since the transformation 6 -^ —9 for all link angles is a symmetry of the Hamiltonian, 
whereas the term sin(6'^i,(x)) changes its sign under this transformation, the expectation value 
(sin(6'^j,(x))) vanishes. Therefore the average energy change of the nonlocal update (|A.3| ) is 



^1 = /5^ E E [1 - coslslV-.+A - ^x))] , (A.5) 

with P = {cos{6fj^i^{x))) = (TrU-p). All expressions for hi that are given in section [^ and in 



section |10| can be obtained from this result by specifying the number of dimensions d and the 
form of the interpolation kernel ip. 

As an example, let us derive eq. ( p.8| ) for abelian gauge theory in two dimensions: For 
piecewise constant interpolation on a two dimensional block we have 

I const J ^ ■'■Ol' ^ fc Xg 



X 



for X ^ Xq 



From eq. ([A.5|) it is clear that only Lb links at the left boundary and Lb links at the right 
boundary of the block contribute to the average energy change hi. This leads to eq. (^.8|): 

hi = 2PPLb[1-cos{s)] . 
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A. 2 Nonabelian gauge fields 

The energy change of a time shce blocking update for SU{2) in d dimensions is 

= -f E E T^^iiKUl.R.+f. - Ul^)Hi:^} , (A.6) 

with if *^^ = U:^+ii^rU*^r,f,Ulr ^nd f/^_^ = gxU:r,f,g*x+f,- H^ is defined analogously. V^ denotes a 
(i — 1 dimensional interpolation kernel, and the rotation matrices R^ G SU{2) are given by 

Rx{n, s) = cos{sipx/'^) + i sin(s-?/'^/2) n-cr . (A. 7) 

Therefore we obtain 

A^ = -f E E I Tr {Ul^m:^) [cos(#./2) cos(#.+^/2) - 1] 

^ xGA[ MT^r *- 



- Tr [in-aUl^Hi:^) sin(#./2) cos(#.+a/2) 
+ Tr (yi^ tn-aHi:^) cos(#./2) sin(#.+^/2) 

- Tr(m-af/f,^2n-<Tii,^;)sin(#,/2)sin(s^,V2)| • (Ai 



Now we will show that to hi = (A7Y) only terms contribute that are even in s. The point is 
that all gauge conditions that are used to define the gf-matrices specify g only up to a global 
gauge transformation: 

gx -^ hgx for all x E Aq . (A. 9) 

If we average the expectation value 

B = (Tr (tn-aUl^H^^:^)) = (Tr {tn-a g^U^^^H^^.g:)) (A.IO) 

over this global symmetry of the Hamiltonian, we get 

B = j dh{Ti{in-ahgxUx,^Hx,^glh*)) = 0. (A.ll) 

Here we used an identity for SU{2) matrices: 

f dhTiiAhBh*) = iTr(A)Tr(5), (A.12) 

and that the Pauli matrices are traceless. 

Therefore the studied updating proposal leads to an average energy change 

hi = PEE I (^Tr(f/.,^ff:J)[l - cos(#./2)cos(#.V2)] 

+ (iTr(2n-af/,^,^^n-aiJf;))sin(s^./2)sin(5^.+^/2)| . (A.13) 
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All expressions for hi that are derived for SU{2) gauge fields in section || and section |T0| can be 
obtained from this formula by specifying the number of dimensions d, the choice of the gauge 
condition that determines the (7-matrices, and the form of the interpolation kernel ip. 

As an example, we derive the expression ( ^.271) for time slice blocking in two dimensions: 
We assume that the ^f-matrices are chosen according to the block axial gauge ( ^.181 ) and that 
i/j vanishes outside the block x'^. Then 



lTr(m-af/,^,^zn-aif,^;) = -iTr (f/,^,^iff 



(A.14) 



since [/|^^ = 1 within the bottom of the block, and we obtain eq. ( p.27| ): 






^x.+a)/2)] • 
Here P = (^Tr([/^^^if^^^)) = (iTr(?7p)) denotes the plaquette expectation value. 
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B Exact results for SU{2) lattice gauge theory in two 
dimensions 



Two dimensional SU{2) lattice gauge theory is defined by the partition function 

Z= fl[dU^,,exp{-niU)). 

•J ^ ,, 



(B.l) 



X,fJ, 



The link variables Ux,^ take values in the gauge group SU{2), and dU denotes the corre- 
sponding invariant Haar measure. The standard Wilson action 'H{U) is given by 



n{U) =/5^[l-iTrf/p] . 

V 

The sum is over all plaquettes in the lattice. Wilson loops are defined by 

W{C) = (iTr(t/(C))) , 



fB.2l 



(B.3) 



where U{C) is the parallel transporter around a Wilson loop C of area A. The exact result for 
W{C) on the infinite lattice is 



W{C) 









(B.4) 



Here, Iv{l3) are the modified Bessel functions. The string tension k is defined by the asymptotic 
behavior of the Wilson loop for large A 



-kA 



W{C) ~ e- 

A^oo 

and has the dimension of a mass squared. Therefore 



The string tension correlation length ^ is given by 

.-1/2 



v^ 



log 



MP). 




^^"U 



(B.5) 



(B.6) 



(B.7) 
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C Gauge conditions on the lattice 

In this appendix we discuss the relation of gauge conditions on the lattice that are formulated 
in terms of parallel transporters Ux,^ with the corresponding gauge conditions given in terms 



of a gauge potential A^^^. We follow the presentation of ref. ^5 . 

Our conventions are as follows: A basis of the SU{N) Lie-algebra is given by traceless, anti- 
hermitean N x N matrices t", normalized to Tr^t^'t'') = — 5"^/2, which satisfy [t"',t''] = /'^''H'^, 
where f"''"^ are the structure constants of the group. In particular for SU{2) we have f^ = i/2(7a, 
where Ua are the Pauli matrices, and /"'"^ = e"'^'^. 

A gauge potential A^^^ on the lattice can be defined by 

^.,. = AIJ'^ = Ip.,, - U:j - ^Tr[f/.,, - U:j . (C.l) 

This definition is motivated by the fact that U^^^ = exp{Ax^fj,) for infinitesimal Aj,^^. 

C.l The Landau gauge 

The Landau gauge on the lattice is defined by p6| 



GL{U,g) =J2J2^eTT {g.U^^.gl^^) ^ maximal . (C.2) 



x£Ko fJ,= l 



By performing the gauge transformation U^^fM — > t/|^^ = gxUx,^ig*c+^ , where g^ is obtained by 
the absolute maximum of Gi{U,g), the absolute maximum is shifted to g^ = 1. The gauge 
configuration U such that Gl{U, g) ioi g = 1 is in an absolute maximum has the property that 
the corresponding gauge potential is transverse: 



d 



a 

^fJ-^x,ii = Z^ [^x+/^,At ~ ^x,ti = • (C-3) 



V^v4"_^ denotes the lattice divergence of A^^^. This can be shown as follows: Consider the one 
parameter subgroup g^r) of the local gauge group defined by 

9{r) = {gx{r)} = {exp(™,)} , oj* = -u . (C.4) 

Here, u^ is an abitrary element of the local Lie algebra of the form u^ = t°-uj^. For fixed u and 
U, let U{t) be the one parameter curve on the gauge orbit through U defined by 

UxA'^) = gx{T)Ux^^gl^^{T) . (C.5) 

Let G{t) be defined by 

G(r) = GL{U,g{r)) = Y. E ^^Tr (exp(™,)f/,,^exp(-™,+^)) . (C.6) 

We have 

dG{T) 



, G'{r) = E 5: Re Tr [(.;,. - a...+^)t/,.,,(r)] 

80 



= Y: EReTrK(f/,,^(r)-f/,._^,^(r))] . 
Now we use the fact that ReTr^UxUx^^ir)) = —ReTr^UxU* Jt)) and get 



G'ir) = 1 5: 5: ReTr [a;.([/.,^(r) - t/;^(r) - Ux-.A^) + U:_,Ar)) 



With the definition (|C.1|) of the gauge potential and taking into account that oOx is traceless, 
we obtain 

G'ir) = E E ReTr [.;.(A.,,(r) - A.-,,,(r))] = -i E ^.% E(^S+,,,(r) - AS,,(r)) . 

The r-dependence of A is induced by the r-dependence of U. Finally, G'{t) can be expressed 
by the lattice divergence ( |C.3|) of the gauge potential: 

G'{r) = -i E ^.%V,A^,^(r) . (C.7) 

xgAo 

If [/ is a stationary point of the gauge functional GL{U,g) at g = 1, we have G'{0) = for all 
cu^, which implies 

V,Al^ = . (C.8) 

Therefore, the transversality of the gauge potential A is the condition for the functional Gl{U, g) 
to be stationary at g = 1. 

C.2 The Coulomb gauge 

The Coulomb gauge condition is 

Gc{U,g)= E ^2^^^"^ {9x^^,^^9*x+^^) = maximal. (C.9) 

Note that now only spatial link variables enter in the gauge functional. A similar argumentation 
as above shows that this gauge condition leads to a transversality condition for the spatial 
components of the gauge potential 



E 

,1=1 



v^l^ = E [^x.+,,. - ^S J = . (c.io) 



Here, V^v4" denotes the d — 1 dimensional lattice divergence of the spatial components A'^ = 
(y4°_^, . . . , A^ ,i_^)^ of the gauge potential. 

C.3 An alternative formulation of the Landau gauge condition 

The Landau gauge condition can also be formulated as a minimization of a quadratic form 

Q{U,g) = {g*,-Ag*) i minimal, (C.ll) 
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under the constraint that the g^ are elements of the gauge group. Here A denotes the gauge 
covariant Laplacian 

( A0). = E [U.,,^^-,, + U:_^J,., - 20 J , (C. 12) 



and the scalar product is defined by 



(0,^)= E ^ReTr(0:^,.)- (C.13) 

The subtraction of a constant of the form 

const = 2 E FReTr((7,(?:) (C.14) 

xsAo 

from Q{U, g) will not alter the solution of the minimization problem. Using this, it is simple to 
check that the gauge conditions ( p.2|) and ( |C.11|) are equivalent. 
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